Two-photon Production in Low-velocity Shocks

Supernova remnants and stellar wind bubbles are iconic examples of shocks in the interstellar medium (ISM). These shocks, with the passage of time, descend to lower velocities. Our interest here is shocks with velocities less than 70 km s⁻¹. The post-shock temperature depends on the mean molecular mass, but we adopt a fiducial value of T_s ≲ 10⁵ K and investigate the cooling of such shock-heated hydrogen gas. These shocks cool primarily via Lyα (whose photons are trapped within the shocked region and eventually die on a dust particle) and two-photon continuum. The latter can be detected by Far Ultra-Violet (FUV) imagers. Low-velocity shocks also arise on Galactic length scales: intermediate-velocity and high-velocity clouds raining down from the lower halo into the disk and gas that is shocked as it enters a spiral arm. Vallée (2017) provides a good description of the Milky Way's spiral arms, and Kim et al. (2008) discuss Galactic interstellar shocks.pre

Stellar bow shocks are another major source of low-velocity shocks. For instance, consider our own Sun, a generic G5V star with a weak stellar wind (2 × 10⁻¹⁴ M_⊙ yr⁻¹ ) moving into a warm (∼7000 K) and partially ionized cloud (ionization fraction, x ≈ 1/3 ) at a relative speed of 23–26 km s⁻¹ (Frisch et al. 2011; McComas et al. 2012; Zank et al. 2013; Gry & Jenkins 2014). Because this velocity may not be larger than the fast-mode magnetosonic velocity of the interstellar cloud, there may be a "bow wave" instead of a bow shock (McComas et al. 2012). In the Galactic disk, interstellar space is occupied by the Warm Neutral Medium (WNM; 10³–8 × 10³ K), the Warm Ionized Medium (WIM; 8 × 10³ K), and the Hot Ionized Medium (HIM; 10⁵–10⁶ K), in roughly equal proportions.

From studies with SDSS-Apogee and Gaia-DR2 (Anguiano et al. 2020), the 3D velocity dispersion of the typical (α-abundance tagged) thin-disk star is 48 km s⁻¹, whereas those belonging to the thick disk have dispersion of 87 km s⁻¹. The majority of these local stars reside in the thin disk with a density ratio n_thin/n_thick = 2.1 ± 0.2. As discussed in a previous study (Shull & Kulkarni 2023), a sizable number of stars should be moving supersonically through ambient gas in the WNM and WIM. ⁴ The sizes of the resulting bow shocks will be determined by the stellar velocity and the magnitude of the stellar wind.

Separately, recent developments warrant a closer look at low-velocity shocks. We draw attention to the discoveries of three large-diameter supernova remnants (Fesen et al. 2021) and a 30° long, thin arc in Ursa Major (Bracco et al. 2020). In large part, these findings were made possible with a new diagnostic—GALEX FUV continuum imaging. The detection of such faint, extended features demonstrates simultaneously the value of the dark FUV sky (O'Connell 1987) as well as the value of the FUV band in detecting two-photon emission, a distinct diagnostic of warm (T ≲ 10⁵ K) shocked gas (Bracco et al. 2020; Kulkarni 2022).

The principal goal of this paper is to develop accurate hydrogen plasma cooling models, paying attention to the production of the two-photon continuum in warm plasma, T ≲ 10⁵ K, the temperature range of interest to low velocity shocks. To this end, we first derive the probability of Lyα, two-photon continuum, and Hα emission resulting from excitation of the ground state of hydrogen to all n ℓ levels for n ≤ 5 (Section 2). Next, we review rate coefficients for line excitation by collisions with electrons (Section 3), followed by a review of collisional ionization (Section 4). The results are combined to construct a cooling curve for warm hydrogen plasma (Section 5). In Section 6 we present a comprehensive (isobaric and isochoric) cooling framework. In Section 7 we apply this model to representative low-velocity shocks and then compute the cumulative photon (Lyα, 2γ and Hα) yields over a grid of post-shock temperatures and initial ionization fractions. In Section 8 we summarize our results and discuss future prospects. Unless otherwise stated, the atomic line data (A-coefficients, term values) are from the NIST Atomic Spectra Database ⁵ and basic formulae from Draine (2011).

Colliding electrons excite hydrogen atoms to various levels and, if sufficiently energetic, ionize H i to H ii. Excited levels are also populated by radiative recombination. Excited hydrogen atoms return to the ground state, some by emitting a Lyman-series photon and others via a cascade of optical/IR recombination lines and ending with Lyα emission. Atoms that find themselves in the metastable 2s²S_1/2 level, if undisturbed over a timescale of ${A}_{2s\to 1s}^{-1}\approx 0.12\,$ s, return to the ground state by emitting a two-photon continuum. Here, A_2s→1s = 8.23 s⁻¹ is the Einstein A-coefficient for the 2s − 1s transition (Drake 1986). Its value should be compared to those for allowed transitions (e.g., 6.26 × 10⁸ s⁻¹ for Lyα and (1–5) × 10⁷ s⁻¹ for Hα, depending on the upper levels, 3s, 3p, 3d, involved).

The goal of this section is to compute the production of Lyα photons, two-photon continuum, and Hα resulting from electronic excitation of H atoms. We consider excitations to 14 n ℓ levels; see Table 2 for term values and index scheme. We make the following two assumptions. (1) The proton density in the plasma is less than the "2s critical density" of 1.5 × 10⁴ cm⁻³ (see Chapter 14 of Draine 2011). This ensures that atoms in the 2s level are not collisionally mixed to the 2p level over a timescale of ${A}_{2s\to 1s}^{-1}$ and thus relax by emitting a two-photon continuum. (2) The cooling plasma is optically thick to Lyman lines, so that Lyman photons are absorbed in the vicinity of where they are emitted. Thus, when computing branching ratios, all allowed Lyman series recombinations can be ignored (case B).

2.1. Photon Yields

Consider, for example, an atom excited to one of the n = 3 levels. An atom excited to 3s or 3d will decay to 2p by emitting Hα followed by Lyα. We ignore forbidden transitions such as ns − 1s two-photon decays; see Chluba & Sunyaev (2008). An atom excited to 3p can decay by emitting Lyβ or decay to 2s by emitting Hα followed by two-photon decay. For the latter, the branching fraction ${{ \mathcal B }}_{\beta }$ for Lyβ emission is A_3p→1s /(A_3p→1s + A_3p→2s ) ≈ 88%. However, under case B, the Lyβ photon will be absorbed elsewhere in the nebula, and the situation will be repeated until de-excitation ends with emission of Hα+Lyα.

For a fiducial value of optical depth (τ_0,α = 1000 ) of Lyα, Table 1 lists the corresponding optical depths for the Lyman series. The branching ratio ${{ \mathcal B }}_{\gamma }$ to emit a Lyγ line is slightly smaller that that for Lyβ. As with Lyβ under case-B conditions, Lyγ will also be converted to some combination of Lyα, optical/IR recombination lines, and a two-photon continuum. The Lyman series oscillator strengths, f ∝ n⁻³, where n is the principal quantum number of the excited state. Thus the Lyman-line optical depths decrease rapidly with increasing n (up the series). In contrast, the branching factors ${ \mathcal B }$ decrease slowly with n.

Table 1. Lyman Lines: Optical Depths and Scatterings

line	λ (Å)	f	τ0	n ℓ	${ \mathcal B }$
Lyα	1215.67	0.4164	1000	2p	1
Lyβ	1025.73	0.07912	160	3p	0.881
Lyγ	972.54	0.02901	56	4p	0.839
Lyδ	949.74	0.01394	26	5p	0.819

Note. Columns 1–4 give the name, wavelength, absorption oscillator strength, and central optical depth of the line. The column density of the nebula is assumed to provide a line-center optical depth of τ₀ = 1000 for Lyα, from which τ₀ for other Lyman lines follow. ${ \mathcal B }$ (column 6) is the branching ratio for an atom excited to an np level (column 5) to relax by emitting the appropriate Lyman series line, as opposed to a multi-decay cascade.

Download table as:  ASCIITypeset image

All states other than ns states have fine structure. For example, the 4p state has two levels, P_1/2 and P_3/2, with very little energy difference between the fine structure levels. The electron collisional excitation rate coefficients presented below (Section 3) refer to the sum of transitions to the entire level, e.g., 1s → 4p. The photon yields for Lyα, 2γ continuum, and Hα are given in Table 2.

Table 2. Photon Yields for Lyα, Hα, and 2γ Continuum

k	Level	pk (Lyα)	pk (Hα)	pk (2γ)
1	2s	0	0	1
2	2p	1	0	0
3	3s	1	1	0
4	3p	0	1	1
5	3d	1	1	0
6	4s	0.585	0.415	0.415
7	4p	0.261	0.261	0.739
8	4d	0.813	0.187	0.187
9	4f	1	1	0
10	5s	0.513	0.378	0.487
11	5p	0.305	0.265	0.695
12	5d	0.687	0.267	0.313
13	5f	0.936	0.702	0.064
14	5g	1	1	0

Note. Photon production yields p_k of Lyα (column 3), Hα (column 4), and two-photon continuum (column 5) upon excitation from ground state (1s) to level n ℓ (column 2). The transitions are also labeled by an index k (column 1). At the precision of this work the energy level depends only on the energy wavenumber: 82,303 cm⁻¹ (n = 2), 97,544 cm⁻¹ (n = 3), 102,879 cm⁻¹ (n = 4) and 105,348 cm⁻¹ (n = 5). As a guide to interpreting this table we note that an H atom excited to k = 3 (3s) relaxes by emitting one Hα photon, one Lyα photon but no two-photon continuum.

Download table as:  ASCIITypeset image

The excitation of lines of hydrogen due to collisions with electrons is a venerable topic in ISM studies. The classic review by Dalgarno & McCray (1972) summarized the atomic physics of the 1960s. Scholz & Walters (1991) undertook detailed calculations of the n = 1 → 2 excitations and also provided an estimate for the H I cooling rate coefficient, Λ_{H I}. Anderson et al. (2002); see also Anderson et al. (2000) present close-coupling R-matrix calculations. We adopt these rates since they offer improved accuracy over previous studies (Scholz et al. 1990). The Anderson et al. (2002) theoretical cross sections were constructed with 15 physical energy states up to n = 5 (1s to 5g) supplemented by 24 pseudo-states described by orbitals ($\overline{n},\overline{{\ell }}$) with $\overline{n}=6\mbox{--}9$ and $\overline{{\ell }}=0\mbox{--}5$.

Anderson et al. (2002) present collision strengths, ${\overline{{\rm{\Omega }}}}_{{ij}}$, for excitation from levels i to j (see Figure 1), averaged over a Maxwellian velocity distribution at electron energies, E_T ≡ k_B T, ranging from 0.5 to 25 eV. The collisional excitation rate coefficients (in cm³ s⁻¹) are then given by:

$$\begin{eqnarray}\begin{array}{rcl}{q}_{i\to j} & = & \displaystyle \frac{2\sqrt{\pi }\ \alpha {a}_{0}^{2}}{{g}_{i}}\sqrt{\displaystyle \frac{{I}_{{\rm{H}}}}{{k}_{B}T}}\,{\overline{{\rm{\Omega }}}}_{{ij}}(T)\,\exp (-{E}_{{ij}}/{kT})\\ & = & \displaystyle \frac{8.629\times {10}^{-6}}{{g}_{i}}\displaystyle \frac{{\overline{{\rm{\Omega }}}}_{{ij}}}{\sqrt{T}}\exp (-{E}_{{ij}}/{k}_{B}T),\end{array}\end{eqnarray} \tag{ 1 }$$

where a₀ is the Bohr radius, α is the fine structure constant, g_i is the degeneracy of level i, and E_ij is the energy difference between level i and j. Since we are interested in excitations from the ground state, we set g_i = 2 for 1s (²S_1/2).

Given our focus on warm plasma, we limit the model fits to 1eV ≤ E_T ≤ 15 eV. After some experimentation, we found that a quadratic polynomial provides an adequate fit ⁶ :

$$\begin{eqnarray}&&{\overline{{\rm{\Omega }}}}_{{ij}}(T)={a}_{0}+{a}_{1}x+{a}_{2}{x}^{2},\end{eqnarray} \tag{ 2 }$$

where $x={\rm{ln}}(T/{10}^{4}\,{\rm{K}})={\rm{ln}}\,{T}_{4}$. The fit is precise to about 1% for all levels except for some 5ℓ levels, for which the fitting errors approach 5%. The fitting coefficients can be found in Table 3.

Table 3. Polynomial Fits to Collision Strengths ^a

k	trans	a0	a1	a2
1	1s-2s	0.2958	0.0074	0.0105
2	1s-2p	0.5422	−0.0819	0.2481
3	1s-3s	0.0688	0.0057	0.0008
4	1s-3p	0.1254	0.0000	0.0382
5	1s-3d	0.0647	0.0052	0.0050
6	1s-4s	0.0246	0.0068	−0.0008
7	1s-4p	0.0466	0.0074	0.0114
8	1s-4d	0.0307	0.0063	0.0016
9	1s-4f	0.0115	−0.0010	0.0002
10	1s-5s	0.0169	0.0022	−0.0004
11	1s-5p	0.0316	−0.0002	0.0055
12	1s-5d	0.0218	0.0016	0.0008
13	1s-5f	0.0090	0.0005	−0.0000
14	1s-5g	0.0042	−0.0012	0.0001

Note.

^a Coefficients (a_i ) for Maxwellian-averaged collision strengths, fitted by ${\bar{{\rm{\Omega }}}}_{k}={\sum }_{i=0}^{2}{a}_{i}{x}^{i}$, where $x=\mathrm{ln}(T/{10}^{4}\,{\rm{K}})$. Here, k is the index of the transition (see Table 2).

Download table as:  ASCIITypeset image

The line cooling rate per unit volume is given by n_e n_{H I}Λ_{H I} where n_e = n_p is the electron (and proton) density and n_{H I} is the density of H atoms. The total particle density is n_t = n_{H I} + n_e + n_p = n_H(1 + x) with n_H = n_p + n_{H I} and x = n_e /n_H.

Zoom In Zoom Out Reset image size

We used the fitting model to compute the run of collisional rate coefficients, q_i→j , with temperature (Figure 2). With the cooling coefficients in hand we computed the sum of the luminosity radiated in lines up to n = 5. We consider this sum to be an adequate representation of Λ_{H I}(T) for warm hydrogen. The H I cooling rate coefficient is

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{\rm{H}}\,{\rm\small{I}}}(T)=\displaystyle \sum _{k=1}^{14}{q}_{k}(T){E}_{k},\end{eqnarray} \tag{ 3 }$$

where E_k is the energy level with index k, the energy of transition with index k (see Table 2). Separately, in Appendix A, we compare this cooling coefficient to previously published coefficients (Spitzer 1978; Scholz & Walters 1991; Dere et al. 1997).

Zoom In Zoom Out Reset image size

The coefficient for energy loss through line X (where, for instance, X denotes Lyα, Hα, 2γ) is given by

$$\begin{eqnarray*}&&{{\rm{\Lambda }}}_{X}(T)=\displaystyle \sum _{k=1}^{14}{q}_{k}(T){p}_{k}(X){E}_{X},\end{eqnarray*}$$

where E_X is line energy and p_k (X) is given in Table 2.

3.1. Simple Fits to Line cooling and Collision Rates

The collisional excitation rate coefficient is the sum over all hydrogen levels,

$$\begin{eqnarray}&&Q(T)=\displaystyle \sum _{k=1}^{14}{q}_{k}(T).\end{eqnarray} \tag{ 4 }$$

Both Q(T) and the hydrogen cooling rate (from excitation to n = 2) fall off with temperature as $\exp (-T/{T}_{12})$, where k_B T₁₂ = 3I_H/4 is the energy difference between n = 1 and n = 2 levels. We fit the collision rate and Λ_{H I} over two temperature ranges: "hot" (10⁴ K < T < 1.5 × 10⁵ K) and "warm" (10⁴ K < T < 1.5 × 10⁴ K ),

$$\begin{eqnarray}&&{Q}_{{\rm{H}}\,{\rm\small{I}}}(T)=A\exp (-T/{T}_{12})\displaystyle \sum _{i=0}^{n}{a}_{i}{x}^{i},\end{eqnarray} \tag{ 5 }$$

where $x=\mathrm{ln}(T/{10}^{4}\,{\rm{K}})$, as before. A similar expression was derived for Λ_{H I}(T). The fitting parameters for Q_{H I} and Λ_{H I} are given in Table 4, and the quality of the fit is displayed in Figure 3.

Zoom In Zoom Out Reset image size

Table 4. Fits to Cooling and Collisional Coefficients

Quantity	A	a0	a1	a2	a3
Λ:hot	6.0 × 10−19	1.018	−0.335	0.290	−0.059
Λ:warm	6.0 × 10−19	1.032	−0.494	0.637
Q:hot	1.0 × 10−7	0.371	−0.132	0.106	−0.021
Q:warm	1.0 × 10−7	0.376	−0.188	0.230

Note. "Quantity" refers to the cooling coefficient (Λ_{H I} in erg cm³ s⁻¹) or collisional excitation rate coefficient (Q_{H I} in cm³ s⁻¹). These quantities are fitted to a model displayed in Equation (5) over two temperature ranges: "hot" (10⁴ K < T < 1.5 × 10⁵ K ) and "warm" (10⁴ K < T < 1.5 × 10⁴ K.)

Download table as:  ASCIITypeset image

We now compute the efficiency of collisional excitations in producing Lyα, Hα, and a 2-photon continuum (each two-photon decay is regarded as one event) as a function of temperature. All inelastic collisional excitations other than to 2s and 2p result in a cascade of photons. The branching ratios are given in Table 2. Both collisional excitation rate coefficients and the branching ratios are used to compute, f_X , the number of photons per collisional excitation; here, X = Lyα, Hα and 2γ. These photon yields are plotted in Figure 4. We see that Lyα has the highest yield, approximately 2/3, followed by 2-photon emission, 1/3. Hα is quite weak, even when measured by photons emitted. This is understandable. Hα emission requires excitation to n ≥ 3 levels, whereas two-photon and Lyα emission are obtained by excitation to n = 2 (and also cascade from higher states). However, Hα has a major advantage: it can be observed with existing ground-based observatories.

Zoom In Zoom Out Reset image size

We also provide model fits for f_X :

$$\begin{eqnarray}&&{f}_{X}=\displaystyle \sum _{k=0}^{m}{a}_{k}{x}^{k},\end{eqnarray} \tag{ 6 }$$

where $x=\mathrm{ln}({T}_{4})$, as before, and m = 2 or m = 1 (as appropriate). The model fit values are displayed in Table 5. The model fit for Hα is displayed in Figure 5.

Zoom In Zoom Out Reset image size

Table 5. Efficiency of Photons per Collision

phot	a0	a1	a2
Hα	0.031	0.131	−0.028
2γ	0.378	−0.042
Lyα	0.622	0.042

Note. Electron collisions with hydrogen atoms produce Lyα, two-photon continuum (2γ), Hα, and other lines. For each of these categories, the efficiency of photon production, f _X , depends on temperature and is fitted to a model given by Equation (6). The fit parameters are displayed above. The model is accurate to 2%; see Figure 5 for an example model; th.

Download table as:  ASCIITypeset image

The collisional ionization rate coefficient is

$$\begin{eqnarray*}&&{k}_{{ci}}(T)={\int }_{{I}_{{\rm{H}}}}^{\infty }{\sigma }_{{ci}}(E)v\,f(E){dE},\end{eqnarray*}$$

where f(E) is the Maxwellian energy distribution, I_H = 13.598 eV is the ionization energy of hydrogen, and σ_ci (E) is the collisional ionization cross section as a function of electron energy in the center-of-mass frame, E = 1/2μ v² with μ = m_e m_H/(m_e + m_H) ≈ m_e . Here, m_e and m_H are the mass of the electron and hydrogen atom, respectively. The collisional ionization rate coefficient can be sensibly written as

$$\begin{eqnarray}&&{k}_{{ci}}(T)=A(T)\exp (-{I}_{{\rm{H}}}/{k}_{B}T).\end{eqnarray} \tag{ 7 }$$

Black (1981) provided an approximate form for k_ci T), based on ionization cross sections tabulated by Lotz (1967),

$$\begin{eqnarray}&&{k}_{{ci}}(T)=5.85\times {10}^{-11}{T}^{1/2}\exp (-{I}_{{\rm{H}}}/{k}_{B}T)\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 8 }$$

This expression is consistent with the approximation, σ_ci ∝ (1 − I_H/E), quoted in Draine (2011), but only valid at low collision energies (I_H ≤ E ≤ 3I_H). As shown by Lotz (1967), the high-energy behavior is ${\sigma }_{{ci}}\propto \mathrm{ln}E/E$. Scholz & Walters (1991) provided a better approximation (for 10⁴ K to 2 × 10⁵ K) using a sixth-order polynomial,

$$\begin{eqnarray}&&A(T)=\exp \left(\displaystyle \sum _{i=0}^{6}{a}_{i}{y}^{i}\right)\ {\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 9 }$$

where $y=\mathrm{ln}T$. A comparison (Figure 6) between Equation (8) (Black 1981 fit) and the more accurate Scholz & Walters fit (Equation (9)) shows that the former breaks down at high temperatures. We offer a modified formula with the correct asymptotic behavior (k_B T > I_H ) as used in the shock models of Shull & McKee (1979),

$$\begin{eqnarray}&&{k}_{{ci}}(T)=\displaystyle \frac{5.85\times {10}^{-11}{T}^{1/2}}{\left[1+0.1{k}_{B}T/{I}_{{\rm{H}}}\right]}\exp \left(-\displaystyle \frac{{I}_{{\rm{H}}}}{{k}_{B}T}\right)\,{{\rm{cm}}}^{3}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 10 }$$

As can be seen from Figure 6, this modified formula provides a good fit at both low and high temperatures. For quick estimates we use Equation (10), but the polynomial formulation of Scholz & Walters (1991; see their Table 2) is preferred when precision is needed (e.g., in numerical integration of differential equations). The resulting ionization power loss per unit volume is n_e n_{H I}Λ_ci where Λ_ci = k_ci I_H is the collisional ionization energy loss coefficient.

Zoom In Zoom Out Reset image size

The goal in this section is to construct a cooling curve for "warm" (T ≲ 10⁵ K) hydrogen plasma. In Section 3 we formulated the cooling coefficient due to line cooling, while in Section 4 we presented the same for ionization losses. Here, we summarize the cooling coefficients for radiative recombination (free-bound) and free–free losses. Armed thus, we formulate the cooling curve for hydrogen in the temperature range 10⁴–10⁵ K.

The kinetic energy of recombining electrons is a loss to the thermal pool. The model fits for radiative effective recombination coefficients α_k where k = 1, A, B (corresponding to recombinations to the n = 1s level, and to all states in case A and case B) can be found in Table 7 (Appendix B). The radiative recombination power loss per unit volume is n_e n_p αΛ_rr where α is either case A or case B, as appropriate. The radiative recombination energy rate coefficient Λ_rr = α〈E_rr〉 where 〈E_rr〉 is the mean thermal energy lost by electron upon recombination. Following Draine (2011) we let 〈E_rr〉 = f_rr k_B T.

The free–free emission rate per unit volume is n_e n_p Λ_ff, where Λ_ff is the free–free emissivity. The free–free power per electron is n_p Λ_ff. The mean time for an electron to recombine is ${({n}_{p}\alpha )}^{-1}$. Thus, the free–free energy lost up to the point of recombination is Λ_ff/α, which we equate to f_ff k_B T. The combined recombination and free–free cooling rate coefficient is then

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{rf}}={{\rm{\Lambda }}}_{\mathrm{rr}}+{{\rm{\Lambda }}}_{\mathrm{ff}}=\alpha {f}_{\mathrm{rf}}(T){k}_{B}T.\end{eqnarray} \tag{ 11 }$$

where f_rf = f_rr + f_ff. The run of f_rf(T) with temperature is displayed in Figure 7 (Appendix B), and the model fits are presented in Table 7 (Appendix B).

Zoom In Zoom Out Reset image size

We now have all the elements to formulate the cooling rate per unit volume, ${ \mathcal C }(T)$, expressed as a negative value (for energy losses):

$$\begin{eqnarray}&&{ \mathcal C }(T)=-{n}_{e}{n}_{{\rm{H}}\,{\rm\small{I}}}\left[{{\rm{\Lambda }}}_{{\rm{H}}\,{\rm\small{I}}}(T)+{k}_{{ci}}(T){I}_{{\rm{H}}}\right]-{n}_{e}{n}_{p}{{\rm{\Lambda }}}_{\mathrm{rf}}.\end{eqnarray} \tag{ 12 }$$

The three RHS terms are given by Equations (3), (7), and (11), respectively.

The investigation of time-dependent cooling of gas heated to T ≈ 10⁵ K is a classic endeavor, constituting the Ph. D. thesis topics of Michael Jura (Jura & Dalgarno 1972) and Minas Kafatos (Kafatos 1973). The motivation in the 1970s seems to have been ambient gas heated by an FUV shock-breakout pulse. Separately, Draine & Salpeter (1978) investigated the production of Lyα from SN shocks, and Shull & Silk (1979) computed UV emission from SNRs in primeval galaxies.

In unmagnetized plasma, assuming a cool (T = 0) pre-shock gas, the post-shock temperature of an adiabatic strong shock (Mach number ≫ 1) is given by

$$\begin{eqnarray}&&\begin{array}{l}{T}_{s}=\displaystyle \frac{2\left(\gamma -1\right)}{{\left(\gamma +1\right)}^{2}}\displaystyle \frac{\mu {v}_{s}^{2}}{{k}_{B}}\\ \,=1.2\times {10}^{5}\,{\rm{K}}\left(\displaystyle \frac{\mu }{{m}_{{\rm{H}}}}\right){\left(\displaystyle \frac{{v}_{s}}{70\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 13 }$$

Here, μ is the mean mass per particle and γ is the ratio of specific heats at constant pressure and constant volume; γ = 5/3 for mono-atomic gas. Let y be the number density of helium relative to that of hydrogen. For the gas at the solar circle, y = 0.1. The mean molecular mass for H⁰ and He⁰ is μ = [(1 + 4y)/(1 + y)]m_H = 1.27m_H. It is μ = 0.67m_H for H⁺ and He⁰, 0.64m_H for H⁺ and He⁺ and 0.61m_H for H⁺ and He⁺².

Three timescales come into play for post-shocked gas: τ_r , the recombination timescale; τ_ci , the collisional ionization timescale; and τ_c , the cooling timescale. For gas around 10⁵ K, we have τ_ci ≪ τ_r . With this inequality, the cooling gas does not obey the conditions for collisional ionization equilibrium. Thus, it is often essential to undertake a full time-dependent calculation.

6.1. Electron-proton Equilibration

At the collisionless shock front, the electrons and protons receive similar amounts of random motion. Being more massive, the protons acquire more energy and are initially much hotter than the electrons. The equilibration timescale for electrons to be heated up to the temperature of the protons via electron-proton encounters is approximately

$$\begin{eqnarray*}&&{t}_{({\rm{e}},{\rm{p}})}=14\ {\left(\displaystyle \frac{T}{{10}^{5}\,{\rm{K}}}\right)}^{3/2}\left(\displaystyle \frac{{\mathrm{cm}}^{-3}}{{n}_{e}}\right)\left(\displaystyle \frac{25}{\mathrm{ln}{\rm{\Lambda }}}\right)\ \mathrm{yr},\end{eqnarray*}$$

where $\mathrm{ln}{\rm{\Lambda }}$ is the Coulomb logarithmic factor accounting for distant encounters (Spitzer 1978; Chapter 2). The current view (Laming et al. 1996; Ghavamian et al. 2007) is that plasma instabilities and electromagnetic waves drive electron-proton equilibration faster than two-body interactions. In Section 6.2 we assume that equipartition occurs before collisional ionization sets in.

6.2. Collisional Ionization

The rate equation for the number density of electrons is

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{e}}{{dt}}={n}_{e}{n}_{{\rm{H}}\,{\rm\small{I}}}{k}_{{ci}}(T)-{n}_{e}{n}_{p}\alpha (T),\end{eqnarray} \tag{ 14 }$$

where α = α_A , α_B as needed. Given our assumption of hydrogen plasma, the number density of protons, n_p = n_e . As noted in Appendix C a solution to this equation at constant T is

$$\begin{eqnarray*}&&x{(t)}^{-1}={x}_{0}^{-1}\exp (-t/{\tau }_{{ci}})+{x}_{{eq}}^{-1}\left[1-\exp (-t/{\tau }_{{ci}})\right].\end{eqnarray*}$$

Here, x₀ = x(t = 0), ${\tau }_{r}\equiv {(\alpha {n}_{{\rm{H}}})}^{-1}$ and ${\tau }_{{ci}}\equiv {({k}_{{ci}}{n}_{{\rm{H}}})}^{-1}$ are the characteristic time scales for recombination and collisional ionization, respectively, and

$$\begin{eqnarray}&&{x}_{{eq}}(T)=\displaystyle \frac{{k}_{{ci}}(T)}{\alpha +{k}_{{ci}}(T)}=\displaystyle \frac{{\tau }_{r}}{{\tau }_{r}+{\tau }_{{ci}}}.\end{eqnarray} \tag{ 15 }$$

The lowest probable value for the ionization fraction in a realistic diffuse atomic medium, before heating commences, is x₀ ≈ 3 × 10⁻⁴. In this case the electrons come from stellar FUV photoionization of C, S, Mg, Si, Fe and other trace metals. As can be seen from Equation (C1) the timescale for electron ionization to reach the half equilibrium fraction is $\sim \mathrm{ln}({x}_{\mathrm{eq}}/{x}_{0}){\tau }_{\mathrm{ci}}$.

6.3. Recombination

As can be seen from Figure 8, collisional ionization is a strong function of temperature. At late times, when the plasma has cooled, collisional ionization can be ignored and Equation (14) simplifies to dx/dt = − x²/τ_r , with the solution

$$\begin{eqnarray*}&&\displaystyle \frac{1}{x}-\displaystyle \frac{1}{{x}_{0}}=\displaystyle \frac{t}{{\tau }_{r}}.\end{eqnarray*}$$

The ionization fraction decreases from x₀ to x₀/m on a timescale of t = (m − 1)τ_r /x₀. The run of recombination timescale as a function of temperature is shown in Figure 8.

Zoom In Zoom Out Reset image size

6.4. Basic Cooling and Recombining Framework

The path in the phase diagram of density, ionization fraction, and temperature along which the gas cools depends on the circumstances. For planar radiative shocks, the pressure behind the shock is ${P}_{0}+(3/4){\rho }_{0}{v}_{0}^{2}$, rising to ${P}_{0}+{\rho }_{0}{v}_{0}^{2}$ downstream when ρ ≫ ρ₀. Here, the pre-shock parameters have subscript 0. Thus, radiative shocks are good examples of cooling at nearly constant pressure ("isobaric"). As the gas cools, its density rises to maintain the pressure. A second possibility is cooling at constant density ("isochoric"). The decrease in temperature, following cooling, leads to lower pressure. Pressure changes are conveyed at the speed of sound, c_s . The timescale for adiabatic sound waves to cross a nebula of length L is τ_a = L/c_s . Isochoric cooling will take place when the cooling time is short, τ_c ≪ τ_a .

The first law of thermodynamics states that any gain in the internal energy (U) of the system is due to increase in internal heat and work done: dU = dQ − PdV. For mono-atomic gas, the internal energy of the nebula per unit volume is U = (3/2)nk_B T, while the pressure is given by Boyle's law P = nk_B T. Here, N = nV is the total number of particles in the nebula whose volume is V. Let N_H = n_H V be the total number of hydrogen nuclei. Ionization can produce changes in N, whereas N_H is fixed.

The three physical parameters governing the cooling hydrogen plasma are n_e , T and n_H. We have two differential equations, one for ionization balance (n_e ; Equation (14)) and one for for energy loss (T; discussed below). A third differential equation follows from the assumed framework: dP/dt = 0 (isobaric cooling) or dV/dt = 0 (isochoric cooling).

For an isochoric system, no work is done by or on the nebula. Adopting case B framework the energy balance equation becomes

$$\begin{eqnarray*}&&q\displaystyle \frac{d}{{dt}}({{nk}}_{B}T)=-{n}_{e}{n}_{{\rm{H}}\,{\rm\small{I}}}\left[{{\rm{\Lambda }}}_{{\rm{H}}\,{\rm\small{I}}}+{k}_{{ci}}{I}_{{\rm{H}}}\right]-{n}_{e}^{2}{\alpha }_{B}{f}_{\mathrm{rf}}{k}_{B}T,\end{eqnarray*}$$

where q = 3/2 and the RHS gives the total cooling rate. Note that n_H remains constant, whereas n = n_H + n_e varies as the ionization fraction changes. The ionization-recombination equation (Equation (14)) can be restated

$$\begin{eqnarray*}&&\displaystyle \frac{{dn}}{{dt}}={n}_{e}{n}_{{\rm{H}}\,{\rm\small{I}}}{k}_{{ci}}-{n}_{e}^{2}{\alpha }_{B}.\end{eqnarray*}$$

We combine the above two equations to obtain

$$\begin{eqnarray*}\begin{array}{rcl}{{qnk}}_{B}\displaystyle \frac{{dT}}{{dt}} & = & -{n}_{e}{n}_{{\rm{H}}\,{\rm\small{I}}}\left[{{\rm{\Lambda }}}_{{\rm{H}}\,{\rm\small{I}}}+{k}_{{ci}}{I}_{{\rm{H}}}+{k}_{{ci}}{{qk}}_{B}T\right]\\ & & -{n}_{e}^{2}{\alpha }_{B}\left({f}_{\mathrm{rf}}-q\right){k}_{B}T,\end{array}\end{eqnarray*}$$

which we deliberately recast as

$$\begin{eqnarray}\begin{array}{rcl}{{qnk}}_{B}\displaystyle \frac{{dT}}{{dt}} & = & -{n}_{e}{n}_{{\rm{H}}\,{\rm\small{I}}}\left[{{\rm{\Lambda }}}_{{\rm{H}}\,{\rm\small{I}}}+{k}_{{ci}}{I}_{{\rm{H}}}\right]-{n}_{e}^{2}{\alpha }_{B}{f}_{\mathrm{rf}}{k}_{B}T\\ & & +{{qk}}_{B}T\left[{n}_{e}^{2}{\alpha }_{B}-{n}_{e}{n}_{{\rm{H}}\,{\rm\small{I}}}{k}_{{ci}}\right].\end{array}\end{eqnarray} \tag{ 16 }$$

In this formulation, the meaning of Equation (16) is clear. The LHS arises from the loss of internal energy. The first term on the RHS represents energy loss from H i collisional line excitation (Λ_{H I}) and collisional ionization (loss of I_H per collision). The loss of kinetic energy per recombination (including the free–free radiation up until the recombination event) is given by the second term. The final term accounts for losses/gains to the thermal pool of electrons during recombination and ionization. For plasma in collisional ionization equilibrium this term vanishes, as expected.

For isobaric cooling, the pressure, P = (n_H + n_e )k_B T is fixed. In this case, we compute n_e and T and then deduce n_H through the pressure equation. As the nebula cools, the ambient gas, in order to maintain the pressure, P_a , does work on the nebula by compressing the nebula. The work done by the medium on the nebula is PdV. However, since P is constant, d(PV) = P_a dV. The internal energy of the nebula is then the enthalpy, HV where H = U + P_a . Going forward, we will drop the subscript to P. It is this store of enthalpy that powers the nebular cooling, ${ \mathcal C }V$. Since (U + P)V = (5/2)Nk_B T we see that Equation (16) still applies but with q = 5/2.

Following a strong shock, the density of the post-shock gas density jumps by a factor of four (for mono-atomic gas, as is the case here). As can be seen from Equation (13) the temperature of the post-shock gas depends on the mean molecular weight of the pre-shocked gas. We elect to specify the initial conditions by physical parameters following the shock: post-shock temperature, T_s and post-shock density H-nuclei density, n_H. The remaining parameter is the ionization fraction. Radiative shock models show that fast shocks, v_s > 110 km s⁻¹, can pre-ionize (H⁺, He⁺) the incoming medium (Raymond 1979; Shull & McKee 1979; Dopita & Sutherland 1996). However, here, we consider shocks with velocities below 100 km s⁻¹. Thus, the ionization fraction of freshly shocked gas is the same as the pre-shocked ionization fraction, x₀.

Consider a representative post-shock temperature, T_s ≈ 10⁵ K. The electron-proton equilibration timescale will be shorter than $14{n}_{e}^{-1}$ yr (see Section 6.1). The collisional ionization time, t_ci, is short, $10{n}_{{\rm{H}}}^{-1}$ yr at T = 10⁵ K, rising to $60{n}_{{\rm{H}}}^{-1}$ yr at 50,000 K. Thus, even if the pre-shock gas has minimal ionization (x₀ ≈ 3 × 10⁻⁴) it will take a time $\sim {\tau }_{\mathrm{ci}}\mathrm{ln}(x/{x}_{\mathrm{eq}})=8{t}_{\mathrm{ci}}$ for the ionization fraction to reach half the equilibrium value. The initial losses are large, owing to both collisional ionization and collisional excitation by the newly liberated electrons and subsequent radiation. Once the gas reaches 10⁴ K, forbidden lines of metals will dominate the cooling. The gas will eventually settle down at T₁ ≈ 5000–8000 K, the temperature of the stable WNM phase (see Heiles & Troland 2003; Kanekar et al. 2003; Murray et al. 2018; Patra et al. 2018). This WNM layer may develop a WIM layer if sufficient diffuse EUV flux is present.

Ignoring "metals," the mean particle mass is μ = m_H(1 + 4y)/(1 + y + x₀). The shock velocity and μ determine the post-shock temperature, T_s through Equation (13). Recall that, in our simplified model, the losses from the shocked nebula are only those associated with hydrogen (line radiation, ionization, free-bound, and free–free). In particular, while we include helium in computing the reduced mass, we do not include losses due to helium. In short, we treat helium as a silent and inactive partner. The energy per H-nucleus and associated electron is E₀ = qk_B T_s (1 + x₀ + y). The end state is when hydrogen has largely recombined and thus the energy per H-atom is E₁ = qk_B T₁.

For our fiducial temperature of 10⁵ K, we have E₀ ≈ [12.9, 21.5](1 + x₀) eV energy per H atom for q = 3/2, 5/2. Thus, on simple grounds, we can see that low-velocity shocks will not significantly affect the ionization of the incoming particles.

7.1. Initial Exploration

We start our investigation by considering two cases, both with post-shock temperature of 1.5 × 10⁵ K and the following parameters: n_H = 100 cm⁻³, x₀ = 3 × 10⁻⁴ (indicative of a shock running into CNM; hereafter case "CNM"); n_H = 1 cm⁻³, x₀ = 0.5 (indicative of the local ISM ("LISM"; Frisch et al. 2011; Linsky et al. 2022). For each case, under the isobaric framework, we compute the evolution of the physical quantities. We plot the resulting run of temperature, electron density and cooling losses with time in Figures 9 and 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to the run of physical quantities (T, x, n_H ). we also plot the cumulative number of recombinations per H nucleus,

$$\begin{eqnarray*}&&{n}_{r}(t)={\int }_{0}^{t}\displaystyle \frac{1}{{n}_{{\rm{H}}}}{n}_{{\rm{H}}}^{2}{x}^{2}{\alpha }_{B}(T){dt},\end{eqnarray*}$$

and the cumulative number of inelastic collisions per atom

$$\begin{eqnarray}&&{n}_{c}(t)={\int }_{0}^{t}\displaystyle \frac{1}{{n}_{{\rm{H}}}}{n}_{{\rm{H}}}^{2}x(1-x)Q(T){dt},\end{eqnarray} \tag{ 17 }$$

where Q(T) is given by Equation (4). A fraction of these inelastic collisions lead to emission of Lyα photons (subscript c=Lyα, n_Ly ), 2-photon continuum emission (n_2γ ), Hα (n_Hα ) photons and ionizations (n_i ).

Radiative recombination also produce Lyα, 2-photon continuum and Hα photons. Draine (2011; Table 14.2) provides, over the temperature range 5 × 10³–2 × 10⁴ K, the recombination coefficients to the 2s level, α_2s and for Hα emission, α_Hα . From these data we find

$$\begin{eqnarray}&&\displaystyle \frac{{\alpha }_{2{\rm{s}}}}{{\alpha }_{B}}\approx 0.328{T}_{4}^{0.115},\ \displaystyle \frac{{\alpha }_{{\rm{H}}\alpha }}{{\alpha }_{B}}\approx 0.450{T}_{4}^{-0.11}.\end{eqnarray} \tag{ 18 }$$

Thus, at typical temperatures of photoionized gas, recombination results in similar diffuse emission for two-photon continuum and Hα. Separately, α_B = α_2s + α_2p. Thus, 1 − α_2s/α_B is the fraction of recombinations which result in production of Lyα photons.

From Figure 9 we see that a shock moving into the CNM (with n_H = 100 cm^-3) cools rapidly. In fact, nominally in less than one year, the post-shocked gas temperature is only 20,000 K. The rapid increase in ionization (panel 2) significantly reduces the proton-ion equilibration time. The strongest energy losses are due to radiative emission (Lyα and 2-photon) followed by loss due to electron collisional ionizations.

Relative to the CNM case, a shock running into the LISM undergoes a slower evolution. At late times (10⁴ yr) as recombination starts the cooling losses mount. This leads to a rapid decrease in temperature and increased recombination. The contribution to the production of line photons (Lyα, 2-photon, Hα) due to electron excitation and recombination is shown in Figure 11. Finally note that, since we are considering low density gas, the photon yields are independent of n_H but do depend on T_s and x₀.

Zoom In Zoom Out Reset image size

7.2. Cumulative Two-photon and Line Yields

A principal goal of the paper is the computation of the number of two-photon decays, Lyα, and Hα per H-nuclei. We distinguish emission due to electron collisions and recombination. The run of these yields as a function of the temperature of the cooling gas is shown in Figure 11. From this figure, we see that the photon "yields" arising from collisions saturate by the time gas has cooled to 10⁴ K, while recombination becomes important only after the temperature has fallen to 2 × 10⁴ K.

These "yields" are a function of shock velocity, v_s and initial ionization fraction, x₀. In Figures 12 and 13 we plot runs of yields due to collisions ("c:N_X " where X = Ly_α , H_α , 2γ ) and due to recombination ("r:N_X ") as a function of x₀ for two representative temperatures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Table 6 we present the calculations for four post-shock temperatures (1.5 × 10⁴ K, 1 × 10⁵ K, 5 × 10⁴ K and 3 × 10⁴ K) as a function of pre-shock ionization fraction, x₀, from 0 to 0.9. Finally, in Figure 14 we display runs of x_f, the ionization fraction by the time gas has cooled to 7 × 10³ K.

Table 6. Photon Yields

Ts	x0	xf	c:NLyα	c:N2γ	c:NHα	r:NLyα	r:N2γ	r:NHα
1.5 × 105	0.00	0.11	1.55	0.70	0.35	0.22	0.03	0.15
1.5 × 105	0.10	0.15	1.70	0.77	0.38	0.29	0.04	0.19
1.5 × 105	0.20	0.18	1.85	0.84	0.41	0.36	0.06	0.25
1.5 × 105	0.30	0.20	2.00	0.90	0.45	0.45	0.07	0.30
1.5 × 105	0.40	0.21	2.13	0.96	0.48	0.53	0.08	0.36
1.5 × 105	0.50	0.21	2.21	0.99	0.50	0.62	0.10	0.41
1.5 × 105	0.60	0.21	2.25	1.01	0.51	0.69	0.12	0.46
1.5 × 105	0.70	0.21	2.27	1.01	0.52	0.75	0.13	0.50
1.5 × 105	0.80	0.21	2.28	1.02	0.52	0.81	0.15	0.53
1.5 × 105	0.90	0.21	2.28	1.02	0.52	0.87	0.16	0.57
1 × 105	0.00	0.05	1.06	0.50	0.22	0.11	0.02	0.08
1 × 105	0.10	0.09	1.16	0.55	0.25	0.17	0.03	0.12
1 × 105	0.20	0.12	1.26	0.59	0.27	0.23	0.04	0.16
1 × 105	0.30	0.15	1.36	0.64	0.29	0.30	0.05	0.20
1 × 105	0.40	0.18	1.46	0.69	0.31	0.36	0.06	0.24
1 × 105	0.50	0.20	1.55	0.73	0.34	0.43	0.07	0.29
1 × 105	0.60	0.21	1.65	0.77	0.36	0.51	0.08	0.34
1 × 105	0.70	0.21	1.71	0.80	0.37	0.58	0.09	0.39
1 × 105	0.80	0.21	1.73	0.81	0.38	0.65	0.11	0.43
1 × 105	0.90	0.21	1.73	0.81	0.38	0.71	0.12	0.47
5 × 104	0.00	0.01	0.53	0.27	0.09	0.03	0.00	0.02
5 × 104	0.10	0.04	0.58	0.30	0.10	0.08	0.01	0.06
5 × 104	0.20	0.07	0.62	0.32	0.11	0.14	0.02	0.09
5 × 104	0.30	0.10	0.67	0.35	0.12	0.19	0.03	0.13
5 × 104	0.40	0.13	0.72	0.37	0.13	0.24	0.04	0.16
5 × 104	0.50	0.15	0.76	0.39	0.14	0.29	0.05	0.20
5 × 104	0.60	0.17	0.81	0.42	0.15	0.35	0.05	0.24
5 × 104	0.70	0.19	0.85	0.44	0.16	0.41	0.06	0.28
5 × 104	0.80	0.20	0.90	0.46	0.17	0.47	0.07	0.32
5 × 104	0.90	0.21	0.93	0.48	0.17	0.54	0.09	0.36
3 × 104	0.00	0.00	0.29	0.16	0.04	0.01	0.00	0.01
3 × 104	0.10	0.03	0.32	0.17	0.04	0.06	0.01	0.04
3 × 104	0.20	0.05	0.34	0.19	0.05	0.11	0.02	0.08
3 × 104	0.30	0.08	0.37	0.20	0.05	0.16	0.02	0.11
3 × 104	0.40	0.11	0.39	0.21	0.06	0.21	0.03	0.14
3 × 104	0.50	0.14	0.42	0.23	0.06	0.26	0.04	0.18
3 × 104	0.60	0.16	0.44	0.24	0.06	0.31	0.05	0.21
3 × 104	0.70	0.18	0.46	0.25	0.07	0.37	0.06	0.25
3 × 104	0.80	0.20	0.48	0.26	0.07	0.43	0.07	0.29
3 × 104	0.90	0.21	0.50	0.27	0.07	0.49	0.08	0.33

Note. The cumulative line and two-photon yields per H-nuclei for post-shock temperature, T_s and initial ionization fraction, x₀ computed in the isobaric framework; see caption to Figure 12 for the symbols used for the yields. The calculations are stopped when gas temperature of 7000 K is reached. x_f is the ionization fraction at the end of the computations.

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

Shocks with velocities near 70 km s⁻¹ abound in our galaxy. Some descend from higher velocity shocks (e.g., supernova remnants) while others start at low velocity (e.g., stellar bow shocks, high velocity cloud shocks). These shocks do not have strong pre-cursor ionization fronts, and as such the post-shocked gas is partially neutral. Such shocks cool primarily through Lyα, two-photon continuum, Hα, and metal emission lines. Lyα is the brightest line, although resonant scattering traps usually traps these photons within the plasma, resulting in absorption by dust grains. Hα is weak but has the great advantage of being observable from the ground.

Two-photon continuum emission is about 50% of Lyα emission (see Table 6). It is several times brighter than Hα, even when one compares photon fluxes rather than energy fluxes. Fortunately, two-photon emission can be observed with space-based observatories. Furthermore, the two-photon continuum has a distinct FUV/NUV ratio. In fact, GALEX FUV and NUV imagery has led to the recent discovery of large middle-aged supernova remnants (Fesen et al. 2021) and exotic shocked stellar bow shocks with angular scales of hundreds of degrees (Bracco et al. 2020). The Ultraviolet Explorer (UVEX) is a NASA Explorer mission currently under development (Kulkarni et al. 2021). Amongst other goals, UVEX aims to undertake FUV and NUV imaging of the entire sky with higher sensitivity and finer spatial resolution, relative to GALEX. The aforementioned successes with GALEX imagery show great promise of identifying and studying low velocity shocks in a future all-sky survey with UVEX.

With this motivation, and using the best available atomic physics data and atomic calculations, we computed the collisional and cooling coefficients for warm hydrogen (T ≲ 10⁵ K). The primary application of the work presented here is in computing two-photon continuum from bow shocks and old supernova remnants. We allow for pre-ionization by keeping the ionization fraction of the pre-shocked gas as a free parameter that can be set to values computed from more sophisticated shock models (Raymond 1979; Shull & McKee 1979; Dopita & Sutherland 1996). Our expectation is that the accurate H-cooling developed here can be incorporated into time-dependent models (e.g., Gnat & Sternberg 2007).

In the Galactic plane and at low latitudes, two-photon emission will be attenuated by dust in the intervening neutral ISM and contaminated by reflected light from dust grains. In practice, this means that the use of two-photon continuum as a diagnostic will be restricted to high Galactic latitudes and will require careful modeling of reflected light. However, the early successes with GALEX (e.g., Bracco et al. 2020; Fesen et al. 2021) promises rich returns from the planned all-sky FUV & NUV survey of UVEX.

We thank Nikolaus Zen Prusinski, California Institute of Technology, for help with CHIANTI, a collaborative project involving George Mason University, the University of Michigan (USA), University of Cambridge (UK) and NASA Goddard Space Flight Center (USA). Partial support was provided to JMS by the astrophysics planning program of New Horizons Extended Mission. We thank the anonymous referee for a most careful reading of our manuscript and excellent suggestions, which resulted in a much better paper.

In this section we compare our cooling curve, Λ_{H I}(T), discussed in Section 5, to three notable published curves.

A.1. Scholz & Walters (1991)

Scholz & Walters (1991) provide collisional ionization coefficients (Equation (9)) and total cooling (line and ionization energy losses) coefficients. The latter is given by the following polynomial model:

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{\rm{H}}}^{{\rm{S}}{\rm{W}}}(T)={10}^{-20}\times \exp \left(\displaystyle \sum _{i=0}^{5}{a}_{i}{y}^{i}-\displaystyle \frac{{T}_{12}}{T}\right){\unicode{x000A0}}{\rm{e}}{\rm{r}}{\rm{g}}\,{{\rm{c}}{\rm{m}}}^{3}\,{{\rm{s}}}^{-1}\end{eqnarray} \tag{ A1 }$$

where T₁₂ = (3/4)T_R with k_B T_R = I_H and $y=\mathrm{ln}T$. We subtracted losses to ionization, Λ_ci = k_ci I_H, to obtain ${{\rm{\Lambda }}}_{{\rm{H}}\,{\rm\small{I}}}^{{\rm{S}}}(T)$. In Figure 15 we compare our Λ_{H I}(T) (Equation (3)) to ${{\rm{\Lambda }}}_{{\rm{H}}\,{\rm\small{I}}}^{{\rm{S}}}(T)$. The Scholz-Walters curve is systematically lower by 15%. After correction for this scale factor, the two curves agree to within ±5%. Separately, Scholz & Walters (1991) codified the collisional coefficients for excitation from the (1s) ground state of hydrogen to f = 2s and 2p levels

$$\begin{eqnarray}&&{q}_{1s\to f}(T)={{\rm{\Gamma }}}_{1s\to f}(T)\exp (-{T}_{12}/T)\end{eqnarray} \tag{ A2 }$$

where Γ_1s→f is formulated as

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}_{1s\to f}(T)=\exp \left(\displaystyle \sum _{i=0}^{5}{b}_{i}{y}^{i}\right){\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}.\end{eqnarray*}$$

In Figure 15 we compare these two coefficients with our own coefficients (Section 3). As with the total line emission cooling, the 2s and 2p collisional rate coefficients are discrepant by the same scale factor.

Zoom In Zoom Out Reset image size

A.2. Spitzer (1978)

A classical formula for the volumetric cooling rate of hydrogen plasma is given in Spitzer (1978), ${ \mathcal C }={n}_{e}{n}_{{\rm{H}}}{{\rm{\Lambda }}}_{\mathrm{Spitzer}}$, where

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{Spitzer}}(T)=7.3\times {10}^{-19}\exp (-{T}_{12}/T)\mathrm{erg}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ A3 }$$

We are aware that our cooling coefficient formulation expresses the line cooling rate per unit volume as n_e n_{H I}Λ_HI. However, there is little difference between our formulation and Spitzer's formulation when the ionization fraction is small. Spitzer's formulation does not include losses from collisional ionization. Spitzer adopted H i cooling rates given in Table 2 of Dalgarno & McCray (1972) which were based on H i excitation cross sections assembled by Gould & Thakur (1970) from theoretical calculations in the 1960s. Those atomic calculations are superseded by more recent work used in this paper (Anderson et al. 2002). In the left panel of Figure 16 we display the ratio of our cooling curves, Λ_{H I} and Λ_H = Λ_HI + k_ci I_H to the cooling rate of Spitzer (Equation (A3)). Our fit is tuned to be accurate over the range 10⁴–10⁵ K. Even bearing this in mind, it appears that Spitzer's formula over-estimates line cooling at temperatures above 10⁴ K by about 25%.

Zoom In Zoom Out Reset image size

A.3. Collisional Ionization Equilibrium: CHIANTI

We conclude this section by briefly discussing hydrogen plasma in collisional ionization equilibrium (CIE); electron ionization balanced by radiative recombination. The volumetric cooling rate is given by ${ \mathcal C }={x}_{\mathrm{eq}}(1-{x}_{\mathrm{eq}}){n}_{{\rm{H}}}^{2}{{\rm{\Lambda }}}_{\mathrm{CIE}}$ where x_eq is given by Equation (15). It is mainly dominated by line cooling. CHIANTI (Dere et al. 1997; Del Zanna et al. 2021) is a major resource for astronomers working on collisionally excited gas, especially hot plasma in CIE. CHIANTI returns x_eq and cooling function ${ \mathcal C }\equiv {n}_{e}{n}_{{\rm{H}}}{\rm{\Gamma }}$ (under assumption of case A). Thus, Γ = (1 − x_eq)Λ_H.

The run of Γ with temperature is displayed in Figure 17. It is evident that Γ from CHIANTI is higher relative to the calculations reported here at higher temperatures. As can be seen from Figure 18 our x_eq agrees with that used in CHIANTI for T > 2 × 10⁴ K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Standard textbooks (e.g., Draine 2011) provide fitting formulae for hydrogen (recombination, free-bound losses) tuned for study of H ii regions. Here, we present fitting formula in the temperature range interest to this paper: 10⁴ K to 10⁵ K. Our starting point is Hummer (1994) who, over an impressive range of 10 K to 10⁷ K, present recombination rate coefficients, α_i where i = 1 (recombination to n = 1) and i = A, B for case A and case B; see Figure 19.

Zoom In Zoom Out Reset image size

Separately, Hummer (1994) also tabulate the kinetic energy loss due to recombination and free–free emission. In Figure 7 we plot the mean kinetic energy, 〈E_rr〉 versus T. Note that 〈E_rr〉 = f_rr k_B T. The resulting fits are presented in Table 7 and are accurate to one percent.

Table 7. Fits to α and f

qty	A	n	b
α1	1.58 × 10−13	−0.518	−0.039
αA	4.16 × 10−13	−0.708	−0.030
αB	2.58 × 10−13	−0.822	−0.045
${f}_{\mathrm{rr}}^{{\rm{A}}}$	0.784	−0.042	−0.020
${f}_{\mathrm{rr}}^{{\rm{B}}}$	0.672	−0.109	−0.021
${f}_{\mathrm{rf}}^{{\rm{A}}}$	1.09	0.035	0.019
${f}_{\mathrm{rf}}^{{\rm{B}}}$	1.17	0.087	0.061

Note. "qty" is fitted to the model of the form $\mathrm{qty}={{AT}}_{4}^{n+b\mathrm{ln}{T}_{4}}$ where the temperature range is 5 × 10³ K < T < 2 × 10⁵ K; here, T₄ = T/(10⁴ K). See text for definition of subscripts. The superscripts stand for case A or case B. The unit for α is cm³ s⁻¹ while that for f is dimensionless. The fits are accurate to 1% over the temperature range of interest.

Download table as:  ASCIITypeset image

For the free–free emission rate coefficient, Λ_ff we use Equation 10.12 of Draine (2011) to derive f_ff(T) (see Section 5 for definition of f_rf ). The values of f_rf = f_rr + f_ff for case A and case B are displayed in Figure 7 (left panel). We assume a simple linear relation and obtain the following fits for f_rf:

$$\begin{eqnarray*}&&\begin{array}{c}{f}_{{\rm{r}}{\rm{f}}}=0.71+0.0154{T}_{4}{\unicode{x000A0}}({\rm{c}}{\rm{a}}{\rm{s}}{\rm{e}}\,{\rm{B}}),\\ {f}_{{\rm{r}}{\rm{f}}}=1.10+0.0898{T}_{4}{\unicode{x000A0}}({\rm{c}}{\rm{a}}{\rm{s}}{\rm{e}}\,{\rm{A}}).\end{array}\end{eqnarray*}$$

We conclude this section with a discussion of collisional ionization. The energy of the electron, E, upon collision goes into ionizing the atom and imparting kinetic energy to the newly liberated electron. We use the low-energy approximation (Equation (8)) for the collisional cross-section, σ_ci ∝ [1 − (I_H/E)] (valid for I_H ≤ E ≤ 3I_H) to find a mean kinetic energy following ionization,

$$\begin{eqnarray*}\begin{array}{rcl}\langle {E}_{k}\rangle & = & \langle E-{I}_{{\rm{H}}}\rangle =\displaystyle \frac{{\displaystyle \int }_{{I}_{{\rm{H}}}}^{\infty }{\sigma }_{{ci}}(E){{vf}}_{E}(E-{I}_{{\rm{H}}}){dE}}{{k}_{{ci}}(T)}\approx 2{k}_{B}T.\end{array}\end{eqnarray*}$$

This formula becomes inaccurate at high temperature (kT ≳ 3I_H) where the adopted fit to σ_ci breaks down (σ_ci peaks and then falls off as $\mathrm{ln}E/E$). The energy 〈E_k 〉 is shared between the colliding electron and the ionized electron. This energy is not a loss since it is returned to the thermal pool. However, over time, the ionized electron will draw qk_B T energy (q = 3/2 in isochoric framework and q = 5/2 in isobaric framework) from the thermal pool. This is a genuine loss (see the discussion following Equation (16)).

Equation (14) can be written as

$$\begin{eqnarray*}&&\displaystyle \frac{{dx}}{{dt}}=\displaystyle \frac{x}{{\tau }_{{ci}}}-\displaystyle \frac{{x}^{2}}{{\tau }_{h}},\end{eqnarray*}$$

where ${\tau }_{h}^{-1}={\tau }_{{ci}}^{-1}+{\tau }_{r}^{-1}$ is twice the harmonic mean of the two timescales. The equilibrium value for ionization is obtained by setting the LHS to zero, x_eq = τ_h /τ_ci . The above equation can be re-arranged to yield

$$\begin{eqnarray*}&&\displaystyle \frac{{{dx}}^{{\prime} }}{{x}^{{\prime} }[1-{x}^{{\prime} }]}=\displaystyle \frac{{dt}}{{\tau }_{{ci}}},\end{eqnarray*}$$

where ${x}^{{\prime} }=x/{x}_{{eq}}$. This equation can be integrated using the method of partial fractions to yield

$$\begin{eqnarray*}&&\mathrm{ln}\left|\displaystyle \frac{{x}^{{\prime} }}{1-{x}^{{\prime} }}\right|=\displaystyle \frac{t}{{\tau }_{{ci}}}+\mathrm{const}.\end{eqnarray*}$$

At t = 0, ${x}^{{\prime} }={x}_{0}/{x}_{{eq}}\equiv {x}_{0}^{{\prime} }$ and thus

$$\begin{eqnarray*}&&\displaystyle \frac{t}{{\tau }_{{ci}}}=\mathrm{ln}\left|\displaystyle \frac{{x}^{{\prime} }(1-{x}_{0}^{{\prime} })}{{x}_{0}^{{\prime} }(1-{x}^{{\prime} })}\right|.\end{eqnarray*}$$

This equation is useful to compute the timescale to achieve a particular level of ionization. As t → ∞ , x → x_eq , as expected. Alternatively, we apply the transformation, u = 1/x:

$$\begin{eqnarray*}&&\left[\displaystyle \frac{{du}}{{dt}}+\displaystyle \frac{u}{{\tau }_{{ci}}}\right]=\displaystyle \frac{1}{{\tau }_{h}}.\end{eqnarray*}$$

Multiply both sides by ${e}^{t/{\tau }_{{ci}}}$ to obtain

$$\begin{eqnarray*}&&\left[{e}^{t/{\tau }_{{ci}}}\displaystyle \frac{{du}}{{dt}}+\displaystyle \frac{u}{{\tau }_{{ci}}}{e}^{t/{\tau }_{{ci}}}\right]=\displaystyle \frac{1}{{\tau }_{h}}{e}^{t/{\tau }_{{ci}}}.\end{eqnarray*}$$

Because the LHS is the derivative of ${{ue}}^{t/{\tau }_{{ci}}}$, the above equation can be readily integrated to yield

$$\begin{eqnarray}&&u(t)={u}_{0}{e}^{-t/{\tau }_{{ci}}}+{u}_{{eq}}\left[1-{e}^{-t/{\tau }_{{ci}}}\right]\end{eqnarray} \tag{ C1 }$$

where $u(t=0)={u}_{0}={x}_{0}^{-1}$ and ${u}_{{eq}}={x}_{{eq}}^{-1}$. As t → ∞ , as expected, u⁻¹ → x_eq .