© ESO 2022

1 Introduction

The 4.8 GHz ortho-formaldehyde (hereafter H₂CO) 1₁₀−1₁₁ masers belong to a small group of rare astrophysical masers. At present only eight Galactic high mass star forming regions with associated 4.8 GHz H₂CO masers are known (Araya et al. 2008; Ginsburg et al. 2015), compared with the more than 900 6.7 GHz methanol masers associated with high mass star forming regions detected in the methanol multibeam survey (Green et al. 2017). Since the H₂CO masers seem to be exclusively associated with high mass star forming regions (Araya et al. 2015), the question arises of what is special about these star forming regions with associated H₂CO masers compared to the majority of high mass star forming regions without associated H₂CO masers.

High luminosity 4.8 GHz megamaser emission has also been detected from the nuclear regions of three ultraluminous infrared (starburst) galaxies (Baan et al. 1986, 2017; Araya et al. 2004, 2008). Analogously to the question about the excitation of the 4.8 GHz masers in a few Galactic high mass star forming regions, the same question arises for the extragalactic H₂CO megamasers. All three knownH₂CO megamaser sources are also OH megamaser sources, which obviously raises the question of whether the two types of masers are excited in the same way.

Since the paper by Boland & de Jong (1981), the question of the pumping of H₂CO masers has regularly been raised over a period of about four decades. Recently, a number of authors (see e.g. Ginsburg et al. 2015; Lu et al. 2019) have remarked that the issue of the pumping of the 4.8 GHz H₂CO masers is still unclear, which makes it difficult to interpret their observations. Further uncertainty also followed from the statement by Baan et al. (2017) that the masers are radiatively pumped, while van der Walt (2014) showed that an inversion of the 1₁₀−1₁₁ transition can be achieved without the presence of an external far-infrared radiation field. van der Walt (2014) also found that no inversion could be found if collisions are switched off.

In view of the uncertainties in the pumping of the H₂CO masers, we expanded the calculations of van der Walt (2014) to better cover the parameter space and thereby to obtain a more complete picture of the inversion of the 1₁₀−1₁₁ transition. We also present the results of calculations aimed at understanding the results obtained with RADEX by Baan et al. (2017) and whether inversion of the 1₁₀−1₁₁ transition can indeed be interpreted as being due to radiative pumping (by a black-body radiation field) as claimed by these authors. Following from the results, we propose a pumping scheme for the H₂CO masers thatexplains why collisions play a central role in the inversion of the 1₁₀ − 1₁₁ transition. For the present purposes we do not consider the effect of beaming, but focus only on the physical conditions under which an inversion can occur.

2 Molecular data

Formaldehyde is a near prolate symmetric top molecule with the dipole moment along the C–O axis (the A-axis). The energy levels are characterized by three quantum numbers, J, K_a, and K_c, with J the total angular momentum, K_a the projection of J on the A-axis (the symmetry axis for a limiting prolate symmetric top), and K_c the projection of J on the C-axis (the symmetry axis for a limiting oblate symmetric top). In the case of o-H₂CO, where the nuclear spins of the two hydrogen atoms are parallel, K_a is odd. The slight asymmetry of the molecule causes each rotational level to be split into two energy levels (known as K-doublets) characterized by different values of K_c.

The energy level diagram for the first 40 levels of o-H₂CO used in this work is shown in Fig. 1 for K_a = 1, 3. The energy levels and Einstein A coefficients for allowed radiative transitions were taken from the Leiden Atomic and Molecular Database (Schöier et al. 2005). The coefficients for collisions of o-H₂CO with o-H₂ and p-H₂, as calculated by Wiesenfeld & Faure (2013), were also obtained from the Leiden Atomic and Molecular Database. In all calculations the ratio of o-H₂ to p-H₂ was that for thermodynamic equilibrium at the specific kinetic temperature. For kinetic temperatures between 100 and 300 K, a linear interpolation was used to produce collision rate coefficients for 10 K intervals instead of the 20 K intervals listed in the Leiden Atomic and Molecular Database.

Formaldehyde has six vibrational normal modes with the lowest at an energy of 1167 cm⁻¹ ≡ 1679 K (Nikitin et al. 2021). This is significantly higher than the excitation energies associated with high mass star forming regions and it is therefore not necessary to consider excitation to the vibrational states.

3 Numerical method

The level populations were found by solving the well-known set of rate equations, which in the escape probability approach is given by $$\begin{array}{lll}\frac{\text{d}{N}_i}{\text{d}t}\hfill & =\hfill & {\sum }_{j<i}[\left(-N_i+\left(\frac{g_i}{g_j}{N}_j-N_i\right)W{\mathcal{N}}_{ij}\right){\beta }_{ij}{A}_{ij}\hfill \\ \hfill & \hfill & +C_{ij}\left(N_j\frac{g_i}{g_j}{e}^{-E_{ij}/kT}-N_i\right)]\hfill \\ \hfill & \hfill & +{\sum }_{j>i}[\left(N_j+\left(N_j-\frac{g_j}{g_i}{N}_i\right)W{\mathcal{N}}_{ji}\right){\beta }_{ji}{A}_{ji}\hfill \\ \hfill & \hfill & +C_{ji}\left(N_j-N_i\frac{g_j}{g_i}{e}^{-E_{ji}/kT}\right)].\hfill \end{array}$$(1)

Here N_i is the number density in level i, g_i the statistical weight of level i, A_ij the Einstein A coefficient for spontaneous emission between levels i and j; W is the geometric dilution factor for an external radiation field with a spectral energy distribution (SED) given by $$F_{\nu }(T_{\text{d}})=[1-e^{-{(\nu /{\nu }_0)}^p}]B_{\nu }(T_{\text{d}}),$$(2)

and ${\mathcal{N}}_{ij}$ the photon occupation number for this field at frequency ν_ij; $C_{ij}=n_{{\text{H}}_2}{K}_{ij}$ is the collision rate with $n_{{\text{H}}_2}$ the H₂ number density and K_ij the collision rate coefficient. The fiducial frequency ν₀ was taken as 3 × 10¹² Hz; β_ij is the escape probability for which we used the expression for the large velocity gradient approximation $${\beta }_{ij}=\frac{1-e^{-{\tau }_{ij}}}{{\tau }_{ij}},$$(3)

with the optical depth given by $${\tau }_{ij}=\frac{A_{ij}}{8\pi }{\left(\frac{c}{\nu }\right)}^3\left(\frac{g_i}{g_j}{x}_j-x_i\right)\frac{N_{\text{col}}}{\Delta v}.$$(4)

The quantities x_i and x_j are the fractional number densities of molecules in the upper level i and lower level j, respectively; N_col is the H₂CO column density and Δv the line width. The ratio N_col∕Δv is the specific column density.

The rate equations were supplemented with the particle number conservation requirement $$N_{\text{tot}}=X{n}_{{\text{H}}_2}={\sum }_i{N}_i,$$(5)

with X the abundance of H₂CO relative to H₂.

We solve the rate equations using a fourth-order Runge-Kutta method. The initial distribution of the level populations was a Boltzmann distribution with temperature equal to the average of the dust and kinetic temperatures. For a given $n_{{\text{H}}_2}$ the calculation started with a small H₂CO specific column density (10⁵ cm⁻³ s) such that the system is completely optically thin in all transitions. An equilibrium solution for this initial specific column density was then found by letting the system evolve in time steps of 5 seconds until the convergence condition |N_i (t_j+1) − N_i(t_j)|∕N_i(t_j) < 10⁻⁶ is reached for all levels. The specific column density is then increased by a small amount (0.001 dex); the equilibrium solution for the previous value of the specific column density is used as the initial distribution. The process is repeated until a specific column density of 5 × 10¹³ cm⁻³ s is reached or is stopped if τ_4.8 > 5 before this specific column density is reached. If there is an inversion in the range of specific column densities, the maximum negative optical depth and the corresponding specific column density are recorded. If there is no inversion τ_4.8 = 0.

Fig. 1 Energy level diagram for the first 40 levels of ortho-formaldehyde. Lower panel: first eight levels with E < 50 K. Upper panel: levels with E > 50 K. The energy scales are different on the y-axes for the lower and upper panels to better separate thedoublet states at lower energies.

4 Results

4.1 The Baan et al. (2017) case

We consider first Fig. 19 in Baan et al. (2017) which shows that inversion of the 1₁₀ − 1₁₁ and 2₁₁ − 2₁₂ transitions of H₂CO can be achieved when locally the gas kinetic temperature is lower than the dust radiative temperature. Understanding the behaviour of the optical depth for different kinetic temperatures is essential before making any conclusions based on the shown behaviour.

Baan et al. (2017) obtained the results shown in their Fig. 19 using the online version of Radex (van der Tak et al. 2007). In particular, it should be noted that the online version of Radex uses an unshielded black-body radiation field (van Langevelde & van der Tak 2008) for the dust radiation. In Fig. 2 we show the dependence of the optical depth on the H₂CO specific column density for the 1₁₀−1₁₁ transition for T_k = 10, 20, 30 K; T_d = 50 K (undiluted black body); and $n_{{\text{H}}_2}$ = 10⁴ cm⁻³ obtained with the above-described numerical method to solve the rate equations. Although there are quantitative differences in the optical depths when compared with Fig. 19 of Baan et al. (2017), the behaviour for the different values of the kinetic temperature is the same, which means that the inversion is largest for T_k = 10 K and becomes smaller as the kinetic temperature increases to approach the dust temperature.

Figure 2 therefore confirms the result of Baan et al. (2017) that under the above-stated conditions the 1₁₀ − 1₁₁ transition shows an inversion. We note that since the radiation field used is that of an undiluted black body, increasing the kinetic temperature to approach the radiation temperature brings the system closer to thermodynamic equilibrium. Thus, when T_k = T_d the level populations follow a Boltzmann distribution, in which case a population inversion is not possible. We illustrate this further in Fig. 3 where we compare the level populations for the lower energy levels where inversions occur, when T_k = 10, 20, 30, 50 K. The inversions for the different transitions are largest for T_k = 10 K and are seen to become smaller as T_k approaches 50 K. For T_k = T_d = 50 K there is no inversion and, consistent with thermodynamic equilibrium, the level populations follow a Boltzmann distribution with a temperature of 50 K.

Still within the framework of the scenario proposed by Baan et al. (2017), we further examined the behaviour of τ_4.8 in the $n_{{\text{H}}_2}$–T_k plane for 10 K ≤ T_k ≤ 300 K, ${10}^4 {\text{cm}}^{-3}\le n_{{\text{H}}_2}\le {10}^6 {\text{cm}}^{-3}$, and when dust emission is that of an undiluted black body at 50 K. The result is shown in Fig. 4. It is seen that there are two distinct regions in the $n_{{\text{H}}_2}$–T_k plane where the 1₁₀−1₁₁ transition is inverted. The first is asmall region with T_k < 40 K and which corresponds with the results shown in Fig. 2. The second is a much larger region with T_k ≳ 130 K and $n_{{\text{H}}_2}$ ≳ 10^4.8 cm⁻³. The absence of an inversion for kinetic temperatures around 50 K is simply a manifestation of the system being close to thermodynamic equilibrium. As we show below, the second region corresponds to the region where inversion is due to collisions only (i.e. without the presence of an external dust infrared radiation field).

In the present case inversion of the 1₁₀−1₁₁, 2₁₁ −2₁₂, and 3₁₂ −3₁₃ transitions has been achieved using a black-body radiation field for the dust emission. In reality the dust radiation field is not that of a black body. Staying within the scenario presented by Baan et al. (2017) that the dust temperature is ~50 K, we used the radio–to–far-infrared SED of Arp 220 derived by Yun & Carilli (2002). The SED of the dust emission is given by Eq. (2) with T_d = 59 K and p = 1.15 (Yun & Carilli 2002) and is compared in Fig. 5 with that of a black body at thesame temperature. In Fig. 6 we show the corresponding variation of |τ_4.8 | in the $n_{{\text{H}}_2}$–T_k plane. Comparison with Fig. 4 shows that while at higher densities and kinetic temperatures the shape of the region where there is an inversion of the 1₁₀ −1₁₁ transition is more or less the same, there is no inversion of the 1₁₀−1₁₁ transition for T_k ≲ 100 K in this case.

Fig. 2 4.8 GHz optical depth as a function of the H2CO specific column density for $n_{{\text{H}}_2}$ = 104   cm−3 and Td = 50 K.

Fig. 3 Comparison of the level populations for levels $J_{K_a,K_c}$ from 11,1 to 41,3 for different kinetic temperatures when $n_{{\text{H}}_2}$ = 104   cm−3 and the external radiation field is a black body with Td = 50 K. The levels are numbered from 1 to 8 just below the case for Td = 50 K.

Fig. 4 Variation in the 4.8 GHz optical depth in the $n_{{\text{H}}_2}$–Tk plane when the spectral energy distribution (SED) of the external emission is given by an undiluted black body at 50 K. τ4.8 = 0 means that there is no inversion.

Fig. 5 Comparison of a black-body SED with the radio–to–far-infrared SED for Arp 220 (Yun & Carilli 2002). The left vertical dashed line is at 140.8 GHz, which is the frequency for the 313 − 111 transition. The right vertical dashed line is at 2808 GHz, which is the highest transition frequency for the 40 levels.

4.2 Expanding the calculations of van der Walt (2014)

van der Walt (2014) presented results of the inversion of the 1₁₀−1₁₁ transition for a limited number of selected combinations of T_k, $n_{{\text{H}}_2}$, and T_d which represented only small parts of the parameter space. We now turn to presenting somewhat more general results covering 10 K ≤ T_k ≤ 300 K and ${10}^4 {\text{cm}}^{-3}\le n_{{\text{H}}_2}\le {10}^6 {\text{cm}}^{-3}$, where we consider the effect of having no external radiation field as a source of excitation and we investigate the effect of an external radiation field on the inversion of the 1₁₀−1₁₁ transition.

In Fig. 7 we show the behaviour of |τ_4.8| in the $n_{{\text{H}}_2}-T_k$ plane when the external radiation field is switched off (i.e. excitation is due only to the internal radiation field and collisions). Inversion of the 1₁₀−1₁₁ transition occurs for T_k ≳ 100 K and $n_{{\text{H}}_2}$ ≳ 10^4.5 cm⁻³. This result supports the results of van der Walt (2014) that the 1₁₀−1₁₁ transition can be inverted without the presence of an external far-infrared radiation field. It does not, however, explain what role collisions play in the pumping of the masers.

With regard to the choice of dust and gas kinetic temperatures, it is noted that Araya et al. (2015) concluded that Galactic H₂CO masers are exclusive tracers of high mass star forming regions. Various studies of Galactic high mass star forming regions have found dust temperatures to be greater than ~50−60 K (see e.g. Kraemer et al. 2001; Barbosa et al. 2016; Lim & De Buizer 2019). This also applies to higher resolution observationsof Arp 220 by Wilson et al. (2014), among others, who derived average dust temperatures of 80 K and 197 K for the eastern and western nuclei, respectively. Downes & Eckart (2007) limits the intrinsic dust temperature in the western nucleus of Arp 220 to between 90 and 180 K based on the size and luminosity of the region. Scoville et al. (2015, 2017) derived a dust temperature >100 K for the western nucleus of Arp 220. We therefore considered three more cases with T_d = 80, 100, and 150 K to illustrate the effects of higher dust temperatures and of geometric dilution. While presenting these results we acknowledge that the combinations of T_d, T_k, and $n_{{\text{H}}_2}$ presented do not necessarily take the thermal balance between dust and gas into account.

In Figs. 8–10 we therefore show the behaviour of |τ_4.8 | in the $n_{{\text{H}}_2}$–T_k plane for T_d = 80, 100, and 150 K, respectively. The lower panels in these figures are for the case of an undiluted radiation field (W = 1 in Eq. (2)). The upper panels in Figs. 8 and 9 are for geometric dilution factors W = 0.3, while the upper panel in Fig. 10 is for W = 0.1. The effect of higher dust temperatures is interesting in that the region in the $n_{{\text{H}}_2}$–T_k plane over which the 1₁₀−1₁₁ transition is inverted changes significantly when W = 1. In Figs. 6 and 7 inversion occurs basically for T_k ≳ 100 K and $n_{{\text{H}}_2}$ ≳   10^4.5 cm⁻³, whereas for warmer dust the inversion occurs at lower H₂ densities and also somewhat lower kinetic temperatures. An undiluted warmer dust radiation field also results in larger values of |τ_4.8| when compared to the case shown in Fig. 7. On the other hand, significant dilution (0.1–0.3) of the external radiation field in all three cases leads to a result very similar to that when collisions and the internal radiation field are the only mechanisms for excitation.

Closer inspection of Figs. 8 and 9 shows the presence of faint regularly spaced vertical regions of weak inversion of the 1₁₀−1₁₁ transition. For example, the top panel of Fig. 9 shows that for W = 0.3 the larger region of inversion connects to lower densities at T_k ~ 140 K and that regions of weak inversions also start to appear at T_k ~ 200 K and 260 K. This results in the wavy nature of the contours for $n_{{\text{H}}_2}$ ≲ 10^4.5 cm⁻³ seen in both Figs. 8 and 9. Test runs with smaller time steps, smaller steps in specific column density, and different convergence criteria were performed to check if this behaviour is due to the numerical method used. In all these cases the wavy behaviour of the optical depth was present. The exact reason for this behaviour of the optical depth is not clear and has not been investigated further.

In addition to the dilution factor, we finally consider the effect of the shape of the SED of the dust emission. In Fig. 11 we show in the upper panel the variation of |τ_4.8| for T_d = 100 K when p = 0.5 in Eq. (2) and in the lower panel when p = 2.0. These two cases should be compared with the lower panel of Fig. 9 for which p = 1.15. It can be seen that although the region over which the 1₁₀−1₁₁ transition is inverted is more or less the same for the three values of p, |τ_4.8 | is smaller over the whole region for p = 0.5 and larger when p = 2.0. In the Rayleigh-Jeans limit and for${(\nu /{\nu }_0)}^p\ll 1$ in Eq. (2), F_ν(T_d) ∝ ν^2+p. When p = 0, F_ν (T_d) = 0.63B_ν(T_d) ∝ ν², which is the flattest possible SED in the Rayleigh-Jeans regime. Except for the factor of 0.63, the result presented in Fig. 4 corresponds to the case when p = 0. Therefore, considering Figs. 4, 9, and 11, there is a clear dependence of |τ_4.8| on the behaviour of the SED at frequencies lower than the turnover frequency: flatter power laws (smaller p) result in smaller values of |τ_4.8| compared to steeper power laws (larger p).

As a last result we note that by following exactly the same calculational procedure as when collisions are included, no evidence could be found for the 1₁₀−1₁₁ transition being inverted by the external dust radiation field (i.e. when collisions are switched off).

Fig. 6 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by $F_{\nu }(T_{\text{d}})=[1-e^{-{(\nu /{\nu }_0)}^p}]B_{\nu }(T_{\text{d}})$ with p = 1.15 and Td = 59 K as for Arp 220 (Yun & Carilli 2002). The synchrotron and free-free components shown in Fig. 5 have also been included.

Fig. 7 Variation of |τ4.8| in the $n_{{\text{H}}_2}-T_{\text{k}}$ plane when there is no external dust radiation field. Excitation is only through collisions and the internal radiation field.

Fig. 8 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by Eq. (2) with p = 1.15 and Td = 80 K. The synchrotron and free-free components shown in Fig. 5 have also been included. The upper panel is for W = 0.3 and the lower panel for W = 1.

Fig. 9 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by Eq. (2) with p = 1.15 and Td = 100 K. The synchrotron and free-free components shown in Fig. 5 have also been included. The upper panel is for W = 0.3 and the lower panel for W = 1.

Fig. 10 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by Eq. (2) with p = 1.15 and Td = 150 K. The synchrotron and free-free components shown in Fig. 5 have also been included. The upper panel is for W = 0.1 and the lower panel for W = 1.

5 Discussion

Inspection of Figs. 7–11 shows some interesting behaviour of the 1₁₀ −1₁₁ maser which had not been seen or anticipated in the results of van der Walt (2014). Although the basic conclusions of van der Walt (2014) are, in a limited sense, still valid, this study shows that the presence of a far-infrared dust radiation field significantly affects the region in the $n_{{\text{H}}_2}$–T_k plane where an inversion can occur. This obviously calls into question the conclusion by van der Walt (2014) that the 4.8 GHz H₂CO masers are strictly collisionally pumped. If a far-infrared radiation field has such a significant effect on the 4.8 GHz masers, then what is the role of collisions, especially if no inversion seems to occur when collisions are switched off? There are two theoretical limiting results which indicate that collisions play a central role in the inversion of the 1₁₀ −1₁₁ maser. The first is that an inversion can be obtained without the presence of an external far-infrared radiation field, and the second that no inversion is found if collisions are switched off even if an external far-infrared radiation field is present. These two limiting cases, as well as the other results, must be the consequence of a specific pumping scheme for the formaldehyde masers. We propose a possible pumping scheme and attempt to explain the results within this framework.

Fig. 11 Comparison of the variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by $F_{\nu }(T_{\text{d}})=[1-e^{-{(\nu /{\nu }_0)}^p}]B_{\nu }(T_{\text{d}})$ with p = 0.5 in the upper panel, p = 2.0 in the lower panel, and Td = 100 K in both cases.

5.1 Proposed pumping scheme

To create a population inversion for the 1₁₀−1₁₁ transition, it is necessary that the excitation rate out of the 1₁₁ state be faster than the rate at which the 1₁₀ state is populated through downward transitions from higher energy states. Inspection of the Einstein A coefficients shows that the fastest radiative downward route to 1₁₀ is through the upper levels of the doublet states of K_a = 1 (e.g. 7₁₆ → 6₁₅ → 5₁₄ → 4₁₃ → 3₁₂ → 2₁₁ → 1₁₀). It therefore seems that the pumping scheme must be such that there is a route and mechanism whereby molecules can be excited to the upper levels of the doublet states. Direct upward radiative excitation out of 1₁₁ to the upper levels of the doublet states is limited to ΔJ = 0, 1,  ΔK_c = ±1 and must therefore start with excitation first to 1₁₀. However, this transition is very slow. For example, for an undiluted dust radiation field as in Eq. (2) with T_d = 100 K and p = 1.15 (the case for Arp 220; the exact value is not important for this example), the radiative transition rate for 1₁₁ → 1₁₀ is only 7.4 × 10⁻¹⁰ s⁻¹. On the other hand, the corresponding rate for 1₁₁ → 2₁₂ is 2.9 × 10⁻⁵ s⁻¹. Upward radiative excitation out of 1₁₁ via 2₁₂ to the lower levels of higher J doublet states (e.g. 3₁₃,   4₁₄,  5₁₅) is therefore much faster than via 1₁₀ to the upper levels of higher J doublet states (e.g. 2₁₁,  3₁₂,  4₁₃,  5₁₄).

To populate the 1₁₀ state, and since radiative excitation out of 1₁₁ via 1₁₀ to the upper levels of the doublet states is very slow, it is necessary to transfer molecules from the lower levels of the doublet states to the upper levels of the doublet states, which then provides a fast radiative decay route to 1₁₀. However, radiative transitions from the lower levels of the doublet states to the associated upper levels are very slow compared to radiative transitions between the lower levels of the doublet states. For example, the radiative transition rate for 5₁₅ → 5₁₄ is 6.2 × 10⁻⁸ s⁻¹ for T_d = 100 K and p = 1.15. On the other hand, the radiative transition rate for 5₁₅ → 4₁₄ is 1.2 × 10⁻³ s⁻¹ and 4.2 × 10⁻⁴ s⁻¹ for 5₁₅ → 6₁₆. The implication is that if a molecule is excited from 1₁₁ to 2₁₂ (the lower doublet state for J = 2), further upward and downward radiative excitations will only be within the ladder of lower level doublet states.

We illustrate the above in Fig. 12, where radiative and collisional rates for different transitions are compared for the case when $n_{{\text{H}}_2}$ = 10^4.25 cm⁻³, T_k = 180 K, T_d = 100 K, and p = 1.15. In the upper panel we show the radiative rates for upward transitions (ΔJ = +1, ΔK_c = +1, black squares) between the lower levels of the doublet states, the corresponding downward rates (ΔJ = −1, ΔK_c = −1, filled red dots), and the upward radiative rates from the lower levels of the doublet states to the corresponding upper levels (ΔJ = 0,  ΔK_c = −1, e.g. 5₁₅ → 5₁₄, open circles) for K_a = 1. Radiative rates between the K_a = 1 and K_a = 3 ladders are so low that radiative transitions between the ladders can be ignored. Also shown are the total collision rates out of the lower level of the doublet states to any other level (blue triangles), and more specifically the collision rates from the lower levels of the doublets to any of the upper level doublet states (red diamonds). As can be seen, the upward radiative rates between the lower levels of the doublet states (black squares) are orders of magnitude higher than the radiative rates from the lower levels of the doublet states to their associated upper levels of the doublet states (open circles). This illustrates the point made above that if a molecule is excited (radiative or collisionally) from 1₁₁ to the 2₁₂ state, further radiative transitions will only be within the ladder of the lower levels of the doublet states and that radiative transfer of molecules from the lower levels of the doublet states to their associated upper levels is not effective. On the other hand, under the given physical conditions, collision rates from the lower levels of the doublet states to any of the upper levels of doublet states (red diamonds) are significantly higher than the corresponding radiative rates. This is also the case, at least for the seven lowest levels of the doublet states above 1₁₁, specifically for collisional transfer of molecules from the lower levels of the doublet states to any of the upper levels of the doublet states (red diamonds). It can therefore be concluded that collisions are the mechanism responsible for population transfer from the lower to the upper levels of the doublet states.

The bottom panel of Fig. 12 shows the corresponding case for the upper levels of the doublet states. Also in this case it is seen that the radiative rates for transitions from the upper levels of the doublet states to the associated lower levels is small compared to radiative rates between the different upper levels of the doublet states and the collisional rates from the upper to the lower levels of the doublet states. Both the upper and lower panels show rather clearly that population exchange between the lower and upper levels of the doublet states is preferentially through collisions and not radiation.

Excitation out of 1₁₁ need not be mainly radiative, but for a given T_k and $n_{{\text{H}}_2}$ it also depends on the energy density in the radiation field at the respective frequencies for the 1₁₁ → 1₁₀ (4.8 GHz) and 1₁₁ → 2₁₂ (~140 GHz) transitions. In Fig. 13 we show the case when p = 2, in other words when the SED for ν ≪ ν₀ (Eq. (2)) is significantly steeper (I_ν ∝ ν⁴) than for p = 1.15 (I_ν ∝ ν^3.15). It is seen that the upward radiative rates for the lower levels of the doublet states (upper panel, black squares) are significantly lower and for the lower J levels even lower than the collisional rates. Collisional excitation out of 1₁₁ actually dominates radiative excitation in this case. The radiative rates for the transfer of molecules from the lower levels of the doublet states to any of the upper levels of the doublet states (open circles) are now also significantly smaller compared to the corresponding collision rates (red diamonds). Thus, as for p = 1.15, population exchange between the lower and upper levels of the doublet states through collisions is by far the dominant process for p = 2.

A possible pumping scheme is therefore that molecules are radiatively or collisionally excited out of the 1₁₁ state. Radiative excitation out of 1₁₁ is almost exclusively to the 2₁₂ state from where further upward radiative excitations to the lower levels of the doublets can take place. Collisional excitation out of 1₁₁ can be either to the lower or upper levels of the doublets. Radiative excitation out of 1₁₁ to 1₁₀ is very slow, which means that the ladder of lower level doublet states will be populated significantly faster than the ladder of upper level doublet states. A fraction of the population in the lower level doublet states is collisionally transferred to the upper level doublet states, which then provides a relatively fast downward radiative route to 1₁₀ to create a population inversion. Figure 14 is a schematic representation of the proposed pumping scheme.

The above-proposed pumping scheme can account, at least qualitatively, for the presented results. First, it is clear that without collisions (i.e. with only the internal and external radiation fields present) there can be no effective population transfer from the ladder of lower doublet states to the upper doublet states to eventually overpopulate the 1₁₀ state relative to the 1₁₁ state. This explains the absence of an inversion when collisions are switched off. Comparison of the bottom panel of Fig. 9 and both panels of Fig. 11 shows that the inversion of the 1₁₀ −1₁₁ transition is smaller for p = 0.5 compared to the cases of p = 1.15 and p = 2. Since the energy density in the radiation field varies as ν^2+p for ν ≪ ν₀, it follows that for p = 0.5 the radiative excitation rates out of the lower levels of the doublets is significantly higher than for p = 2. The upper panel of Fig.12 also shows that the upward and downward radiative rates for the lower levels of the doublet states are very similar for p = 1.15. As a consequence, molecules are more or less ‘trapped’ within the ladder of lower doublet states when p = 0.5 with the probability to be collisionally transferred to the ladder of upper levels of the doublet states being small. On the other hand, for p = 2 the downward rates within the ladder of lower level doublet states are markedly higher than the upward rates, meaning that molecules in higher J states will rather quickly radiatively decay to lower states where the probability for collisional transfer to the ladder of upper level doublet states is comparable to that of upward radiative excitation within the ladder of the lower levels of the doublet states. Such a scenario can explain why the inversion is larger for p = 2 than for p = 0.5.

Fig. 12 Radiative and collision rates of the different energy levels for the lower levels of the doublet states (upper panel) and upper levels of the doublet states (bottom panel). The rates were calculated for Td = 100 K, Tk = 180 K, $n_{{\text{H}}_2}$ = 104.25 cm−3, and p = 1.15. Radiative up (black squares) means upward radiative excitation out of a state (absorption); radiative down (red filled circles) means the sum of spontaneous and stimulated emission rates; collisions (blue triangles) means total collision rate out of a state; radiative lower to upper (open circles) means radiative rate for absorption from the lower level of a doublet state to the associated upper level; collisions lower to upper (red diamonds) means the collision rate from the lower level of a doublet state to its associated upper level.

Fig. 13 Same as for Fig. 12, but for p = 2.0.

Fig. 14 Energy level diagram of lower levels of the doublets and the upper levels of the doublets to illustrate the proposed pumping scheme for the 4.8 GHz H2CO masers. Collisions couple the two sets of states.

5.2 Constraints due to the H₂CO abundance

Figures 8 to 11 suggest that inversion of the 1₁₀ −1₁₁ transition can happen over a significant part of the parameter space with the implication that the 4.8 GHz masers should be quite abundant. However, the H₂CO masers are very rare, which means that there must be other factors that affect the occurrence of these masers. van der Walt (2014) noted that his results imply that the H₂CO abundance must be significantly higher than what is generally observed. In our calculations the H₂CO abundance, $X_{{\text{H}}_2\text{CO}}$, is implicitly included in the specific column density given by $X_{{\text{H}}_2\text{CO}}{n}_{{\text{H}}_2}\mathcal{l}/\Delta v$, where ℓ is the maser path length and Δv the maser line width. In the variation of τ_4.8 with the specific column density (e.g. as in Fig. 2) the most obvious specific column density that can be used to calculate the abundance is that where the optical depth has its maximum negative value. Figure 15 shows, as an example, the variation of the optical depth and $X_{{\text{H}}_2\text{CO}}$ at T_k = 180 K for T_d = 80, 100, 150 K, W = 1 (Figs. 8– 11), Δv = 10⁵ cm s⁻¹, and ℓ = 10¹⁷ cm (Cragg et al. 2002). The variation of the optical depth with $n_{{\text{H}}_2}$ for the three dust temperatures is seen to be quite different. On the other hand, the derived H₂CO abundances behave in the same way and vary between ~1.0 × 10⁻⁴ and ~1.3 × 10⁻⁷ from $n_{{\text{H}}_2}$ = 10⁴   cm⁻³ to 10⁶   cm⁻³. Using the filling factor corrected data presented in Table 8 of Ginsburg et al. (2011), an average $n_{{\text{H}}_2}$ of 4.5 × 10⁴ cm⁻³ and average H₂CO abundance of 1.1 × 10⁻¹¹ for more than 60 lines of sight were found. It is clear that the observed abundances are orders of magnitude too small toexplain the masers. The required abundances are significantly higher than most models can produce (see e.g. Guzmán et al. 2011). It is worth noting, however, that Maret et al. (2004) derived a H₂CO abundance of 6 × 10⁻⁶ for the warm inner envelope of the young low mass star L1527. This suggests that similar cases may also exist for young high mass stars.

Recently Vichietti et al. (2016) proposed a chemical route for the production of H₂CO in young high mass star forming regions leading to abundances as high as 10⁻⁵. Using this as a guide, we illustrate in Fig. 16 what the effect is of setting an upper limit of 10⁻⁵ for the H₂CO abundance on the region in the $n_{{\text{H}}_2}$–T_k plane wherean inversion can occur. Comparison with Fig. 9 shows that the allowed region for an inversion is significantly smaller than when no constraint is put on the H₂CO abundance. Such a significantly smaller region in the $n_{{\text{H}}_2}$–T_k plane might be a reason for the rarity of the H₂CO masers.

An important aspect of Fig. 16 is that for W = 0.3 (upper panel) inversion of the 1₁₀ −1₁₁ transition occurs basically for $n_{{\text{H}}_2}>{10}^{5.5}{\text{cm}}^{-3}$. Within the framework of our calculations the implication is that the collision rates as shown in Figs. 12 and 13 (blue triangles and red diamonds) have to be increased by a factor of 1.25 dex. If this is done the collisional excitation rate out of 1₁₁ is seen to dominate the radiative excitation rate by more than an order of magnitude. If this is the case it can indeed be said that the H₂CO masers are collisionally pumped. In fact, this explains why there is a region in the $n_{{\text{H}}_2}$–T_k plane where an inversion can occur without the presence of an external radiation field.

Fig. 15 Plot showing $X_{{\text{H}}_2\text{CO}}$ (solid lines) and 4.8 GHz optical depth (dashed lines) as a function of H2 density for Tk = 180 K, p = 1.15, and W = 1 for different dust temperatures. A maser path length of 1017 cm and line width of 1 km s−1 was used to calculate the H2CO abundance.

Fig. 16 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane, with the SED of the dust emission given by Eq. (2), with p = 1.15 and Td = 100 K and when the H2CO abundance is restricted to be less than 10−5. The upper panel is for W = 0.3 and the lower panel for W = 1. This result should be compared with that shown in Fig. 9 where no upper limit is placed on the H2CO abundance.

5.3 Evaluation of the results of Baan et al. (2017)

As is shown in Sect. 4.1, using the same set of parameter values as was used by Baan et al. (2017), the population inversion became weaker when T_k approached T_d and disappeared completely when T_k = T_d. As we show in Fig. 3, the disappearance of the population inversion is most definitely due to the system being in thermodynamic equilibrium when T_k = T_d and that a population inversion is not possible in such a case. It can also be seen, especially from Figs. 9 and 10, that the condition T_k < T_d as a requirement for an inversion, as set by Baan et al. (2017) (and not by van der Walt 2014, as stated by these authors), does not follow from our results.

The results presented above have some definite implications for interpreting the H₂CO megamaser emission presented by Baan et al. (2017). It is first necessary to note that no inversion is found for T_k ≲ 100 K when using a more realistic SED with T_d = 59 K for the dust emission (Fig. 6). A direct implication of this result is that, contrary to the conclusions by Baan et al. (2017), the 4.8 GHz H₂CO megamasers in the three starburst galaxies are not associated with cold material. Furthermore, if the effect of an upper limit of 10⁻⁵ on the abundance of o-H₂CO is taken intoaccount (Fig. 16), the implication is that the masers trace hot (100–300 K) gas with densities 10⁵ −10⁶ cm⁻³ rather than cold gas. Such kinetic temperatures have recently been reported by Gieser et al. (2021) for quite a number of high mass star forming regions. It is fairly well established that the 4.8 GHz masers in the Galaxy are associated exclusively with high mass star forming regions (Araya et al. 2015). There is thus no reason why it should not also be the case for starburst galaxies. In the three starburst galaxies observed by Baan et al. (2017) the 4.8 GHzmaser emission is spatially extended and associated with star forming regions, which by itself suggests the association of these masers with high mass star forming regions as in the case of the Galaxy. Why H₂CO maser emission is associated with only three starburst galaxies, and exactly why it is rare in the Galaxy, are questions still to be answered.

6 Conclusions

We performed a new set of numerical calculations to investigate the inversion of the 1₁₀ −1₁₁ transition of o-H₂CO. The calculations covered a significantly larger part of parameter space compared to previous calculations. Considering the results presented above we summarize and conclude as follows:

• 1.

  While theoretically there is a population inversion of the 1₁₀−1₁₁ transition when using a 50 K black-body radiation field, the decrease in the optical depth as the kinetic temperature is increased from 10 to 40 K, as presented if Fig. 19 of Baan et al. (2017), is due to the system approaching thermodynamic equilibrium. This behaviour does not imply that the inversion is the consequence of far-infrared pumping. When a more realistic grey-body SED corresponding to that of Arp 220 is used, no inversion of the 1₁₀ −1₁₁ transition is found for T_k < 100 K. Our calculations suggest that the megamasers in the three starburst galaxies observed by Baan et al. (2017) are not associated with cold material, but rather with dense and hot gas associated with high mass star formation.

• 2.

  We show that, within the framework of our calculations, the 1₁₀−1₁₁ transition can be inverted without the presence of a far-infrared radiation field. On the other hand, no evidence could be found for the 1₁₀−1₁₁ transition being inverted by the external dust radiation field (i.e. when collisions are switched off). Collisions therefore seems to play a key role in the pumping mechanism of the masers.

• 3.

  We proposed a pumping scheme in which molecules are excited either radiatively or collisionally out of the 1₁₁ state. Radiative excitation out of 1₁₁ leads to the population of the ladder of lower levels of the doublet states. Population exchange between the ladders of lower and upper doublet states is collisional since radiative coupling between the two ladders is weak. At high H₂ densities collisional excitation out of 1₁₁ can be to both ladders. Molecules excited collisionally to higher states in the ladder of lower levels of the doublet states still have to be transferred collisionally to the ladder of upper levels of the doublet states. Without collisions, population exchange between the ladders of the lower and upper levels of the doubletstates is extremely slow. Transfer of molecules from the lower levels of the doublet states to the upper levels of the doublet states is necessary since it then provides a fast radiative decay route to 1₁₀.

• 4.

  When assuming reasonable values for the maser line width and the maser path length it is found that, over most of parameter space where theoretically an inversion can occur, the required abundance of H₂CO is significantly larger than observed and predicted by models. Limiting the H₂CO abundance to less than 10⁻⁵ significantly reduces the region in the $n_{{\text{H}}_2}$–T_k plane where an inversion can occur to a rather narrow region with H₂ densities greater than ~10⁵ cm⁻³ and T_k ≳ 100 K. This might explain the rarity of the H₂CO masers in the Galaxy and starburst galaxies.

We finally pointout that o-H₂CO has several submillimeter transitions between excited states (for which E_upper > 100 K) which can be used to probe the physical conditions of the environments of the H₂CO masers. Dedicated high resolution observations of these transitions can be used to test the predictions of our model calculations.

Acknowledgements

This work is based on research supported in part by the National Research Foundation of South Africa (Grant Numbers 132494).

References

All Figures

	Fig. 1 Energy level diagram for the first 40 levels of ortho-formaldehyde. Lower panel: first eight levels with E < 50 K. Upper panel: levels with E > 50 K. The energy scales are different on the y-axes for the lower and upper panels to better separate thedoublet states at lower energies.
In the text

	Fig. 2 4.8 GHz optical depth as a function of the H2CO specific column density for $n_{{\text{H}}_2}$ = 104   cm−3 and Td = 50 K.
In the text

	Fig. 3 Comparison of the level populations for levels $J_{K_a,K_c}$ from 11,1 to 41,3 for different kinetic temperatures when $n_{{\text{H}}_2}$ = 104   cm−3 and the external radiation field is a black body with Td = 50 K. The levels are numbered from 1 to 8 just below the case for Td = 50 K.
In the text

	Fig. 4 Variation in the 4.8 GHz optical depth in the $n_{{\text{H}}_2}$–Tk plane when the spectral energy distribution (SED) of the external emission is given by an undiluted black body at 50 K. τ4.8 = 0 means that there is no inversion.
In the text

	Fig. 5 Comparison of a black-body SED with the radio–to–far-infrared SED for Arp 220 (Yun & Carilli 2002). The left vertical dashed line is at 140.8 GHz, which is the frequency for the 313 − 111 transition. The right vertical dashed line is at 2808 GHz, which is the highest transition frequency for the 40 levels.
In the text

	Fig. 6 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by $F_{\nu }(T_{\text{d}})=[1-e^{-{(\nu /{\nu }_0)}^p}]B_{\nu }(T_{\text{d}})$ with p = 1.15 and Td = 59 K as for Arp 220 (Yun & Carilli 2002). The synchrotron and free-free components shown in Fig. 5 have also been included.
In the text

	Fig. 7 Variation of |τ4.8| in the $n_{{\text{H}}_2}-T_{\text{k}}$ plane when there is no external dust radiation field. Excitation is only through collisions and the internal radiation field.
In the text

	Fig. 8 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by Eq. (2) with p = 1.15 and Td = 80 K. The synchrotron and free-free components shown in Fig. 5 have also been included. The upper panel is for W = 0.3 and the lower panel for W = 1.
In the text

	Fig. 9 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by Eq. (2) with p = 1.15 and Td = 100 K. The synchrotron and free-free components shown in Fig. 5 have also been included. The upper panel is for W = 0.3 and the lower panel for W = 1.
In the text

	Fig. 10 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by Eq. (2) with p = 1.15 and Td = 150 K. The synchrotron and free-free components shown in Fig. 5 have also been included. The upper panel is for W = 0.1 and the lower panel for W = 1.
In the text

	Fig. 11 Comparison of the variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane when the SED of the dust emission is given by $F_{\nu }(T_{\text{d}})=[1-e^{-{(\nu /{\nu }_0)}^p}]B_{\nu }(T_{\text{d}})$ with p = 0.5 in the upper panel, p = 2.0 in the lower panel, and Td = 100 K in both cases.
In the text

Fig. 12 Radiative and collision rates of the different energy levels for the lower levels of the doublet states (upper panel) and upper levels of the doublet states (bottom panel). The rates were calculated for Td = 100 K, Tk = 180 K, $n_{{\text{H}}_2}$ = 104.25 cm−3, and p = 1.15. Radiative up (black squares) means upward radiative excitation out of a state (absorption); radiative down (red filled circles) means the sum of spontaneous and stimulated emission rates; collisions (blue triangles) means total collision rate out of a state; radiative lower to upper (open circles) means radiative rate for absorption from the lower level of a doublet state to the associated upper level; collisions lower to upper (red diamonds) means the collision rate from the lower level of a doublet state to its associated upper level.

In the text

	Fig. 13 Same as for Fig. 12, but for p = 2.0.
In the text

	Fig. 14 Energy level diagram of lower levels of the doublets and the upper levels of the doublets to illustrate the proposed pumping scheme for the 4.8 GHz H2CO masers. Collisions couple the two sets of states.
In the text

	Fig. 15 Plot showing $X_{{\text{H}}_2\text{CO}}$ (solid lines) and 4.8 GHz optical depth (dashed lines) as a function of H2 density for Tk = 180 K, p = 1.15, and W = 1 for different dust temperatures. A maser path length of 1017 cm and line width of 1 km s−1 was used to calculate the H2CO abundance.
In the text

	Fig. 16 Variation of |τ4.8| in the $n_{{\text{H}}_2}$–Tk plane, with the SED of the dust emission given by Eq. (2), with p = 1.15 and Td = 100 K and when the H2CO abundance is restricted to be less than 10−5. The upper panel is for W = 0.3 and the lower panel for W = 1. This result should be compared with that shown in Fig. 9 where no upper limit is placed on the H2CO abundance.
In the text