NEBULAR: A Simple Synthesis Code for the Hydrogen and Helium Nebular Spectrum

Free-bound emission is caused by ions capturing free electrons. Ercolano & Storey (2006) have tabulated the continuous free-bound recombination emission coefficients, γ_bf, of H i, He i and He ii (the ionic states after recombination) over temperatures 10² ≤ T ≤ 10⁵ K with a step size of Δ log₁₀T = 0.1 (hereafter, the logarithmic bases are omitted). NEBULAR performs a cubic spline interpolation over temperatures, and a linear interpolation over frequencies within discontinuities. As stated in Ercolano & Storey (2006), bilinear interpolation reproduces the true values to within 1% or better.

Free–free emission is caused by electrons scattering in the Coulomb field of ions with nuclear charge Z. NEBULAR computes the free–free spectrum for H⁺, He⁺ and He⁺⁺. The emission coefficients are

$$\begin{eqnarray}{j}_{\mathrm{ff}}(\nu ) & = & \displaystyle \frac{1}{4\pi }\,{N}_{X}\,{N}_{e}\ \displaystyle \frac{32\,{Z}^{2}{e}^{4}h}{3{m}_{e}^{2}{c}^{3}}{\left(\displaystyle \frac{\pi h{\nu }_{0}}{3{k}_{{\rm{B}}}T}\right)}^{1/2}\\ & & \times \,\exp (-h\nu /{k}_{{\rm{B}}}T)\ {g}_{\mathrm{ff}}(Z,T,\nu )\end{eqnarray} \tag{ 5 }$$

(e.g., Brown & Mathews 1970). Here, e, m_e, h, and k_B are the electron charge, electron mass, Planck constant and Boltzmann constant, respectively. hν₀ denotes the ionization energy for hydrogen. Note that this equation must be evaluated in cgs units, in particular the electronic charge is 1e = 4.803207 × 10⁻¹⁰ statC, with 1 statC = 1 cm^3/2 g^1/2 s⁻¹.

g_ff(Z, T, ν) is the Gaunt factor, a quantum mechanical correction to the classical scattering formalism. For typical conditions in astrophysical nebulae, g_ff ∼ 1−2. Non-relativistic Gaunt factors have been calculated by van Hoof et al. (2014) for a wide range of temperatures, frequencies and nuclear charges. NEBULAR uses the tabulated temperature-averaged $\langle {g}_{\mathrm{ff}}({\gamma }^{2},u)\rangle$ for Maxwell–Boltzmann electron distributions, where ${\gamma }^{2}={Z}^{2}\,h{\nu }_{0}/({k}_{{\rm{B}}}T)$ and u = hν/(k_BT).

The data for $\langle {g}_{\mathrm{ff}}({\gamma }^{2},u)\rangle$ are available on a log-log grid with regular step sizes in both dimensions. NEBULAR interpolates this grid using a 2D bicubic interpolation from the GNU Scientific Library (GSL v2.1).

2.3.1. H i and He ii

The 2²S level in hydrogenic ions can e.g., be populated by direct recombination onto that level, and by downward cascades following electron capture to a higher level. A radiative decay to the 1²S level is forbidden by the dipole selection rules (Laporte & Meggers 1925). However, the transition may occur if two photons are emitted whose energies add up to the Lyα energy. The transition probability for the hydrogenic $2{\rm{S}}\to 1{\rm{S}}$ two-photon decay is

$$\begin{eqnarray}&&{A}_{2q}=8.2249\,{Z}^{6}\,\displaystyle \frac{{R}_{Z}}{{R}_{H}}\,{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 6 }$$

(Nussbaumer & Schmutz 1984). Here, R_H and R_Z are the Rydberg constants for hydrogen and for hydrogenic ions with nuclear charge Z, respectively. Two-photon decays from higher levels (Hirata 2008) are ignored, as they contribute on the sub-percent level, only.

The two-photon spectrum has a natural upper cut-off at the Lyα frequency. It peaks at half the Lyα frequency when expressed in frequency units, emphasizing the importance of this process for restframe near-UV and far-UV observations. For low densities, the hydrogen two-photon emission dominates the hydrogen free-bound continuum (Section 2.1) below 2000 Å.

The two-photon emission coefficient is

$$\begin{eqnarray}&&{\gamma }_{2q}(\nu )={\alpha }_{{2}^{2}{\rm{S}}}^{\mathrm{eff}}(T,{N}_{e})g(\nu ){\left(1+\displaystyle \frac{{\displaystyle \sum }_{X}{N}_{X}{q}^{X}}{{A}_{2q}}\right)}^{-1}.\end{eqnarray} \tag{ 7 }$$

Here, ${\alpha }_{{2}^{2}{\rm{S}}}^{\mathrm{eff}}$ is the effective recombination coefficient, g(ν) the frequency dependence, and q^X the collision rate transition coefficient; see also Equation (4.29) in Osterbrock & Ferland (2006), and Equation (9) in Storey & Sochi (2014).

The three factors are discussed in the following.

2.3.1.1. Effective Recombination Coefficient ${\alpha }_{{2}^{2}{\rm{S}}}^{\mathrm{eff}}$

${\alpha }_{{2}^{2}{\rm{S}}}^{\mathrm{eff}}(T,{N}_{e})$ is the effective recombination coefficient to the 2²S level. It has been tabulated for both H i and He ii by Hummer & Storey (1987) over wide temperature and density ranges (Case B, only). The coefficients from the original publication were extracted and digitized, and are reproduced in Appendix C (Tables 2 and 3). For convenience, these tables contain the coefficients for the 2²P level as well, so that the total n = 2 effective recombination coefficient may be calculated (not needed for the two-photon emission). For H i, NEBULAR extrapolates ${\alpha }_{{2}^{2}{\rm{S}}}^{\mathrm{eff}}(T,{N}_{e})$ from log T = 3.0 to 2.5, from log T = 4.5 to 5.0, and from log N_e = 2.0 to 0.0. For He ii, extrapolation is used from log T = 3.5 to 3.0, and from log N_e = 2.0 to 0.0. These extrapolations are expected to be accurate within a few percent, as the functional dependence of ${\alpha }_{{2}^{2}{\rm{S}}}^{\mathrm{eff}}(T,{N}_{e})$ at these densities and temperatures is weak.

2.3.1.2. Frequency Dependence g(ν)

The frequency dependence of the two-photon spectrum is absorbed in

$$\begin{eqnarray}&&g(\nu )=h\,\displaystyle \frac{\nu }{{\nu }_{12}}\,A\left(\displaystyle \frac{\nu }{{\nu }_{12}}\right)/{A}_{2q}.\end{eqnarray} \tag{ 8 }$$

Here, ν₁₂ is the Lyα frequency of the respective hydrogenic ion. Nussbaumer & Schmutz (1984) have presented an analytic fit for the frequency dependent transition probability, A(ν/ν₁₂). It is accurate to within 1% over nearly the full spectral range,

$$\begin{eqnarray}&&A(\nu /{\nu }_{12})=C[y(1-{(4y)}^{\gamma })+\delta {y}^{\beta }{(4y)}^{\gamma }],\end{eqnarray} \tag{ 9 }$$

with $y=\tfrac{\nu }{{\nu }_{12}}\left(1-\tfrac{\nu }{{\nu }_{12}}\right)$, β = 1.53, γ = 0.8, δ = 0.88 and C = 202.0 Z⁶ R_Z/R_H s⁻¹, where Z is the nuclear charge.

2.3.1.3. Collisional De-population of the 2²S Level

The 2²S level can be de-populated by means of angular momentum changing (l-changing) collisions, ${2}^{2}{\rm{S}}\to {2}^{2}{\rm{P}}$, inflicted by other ions and electrons. The 2²P level may then decay to the 1²S level by emitting a Lyα photon, reducing the two-photon flux. This is described by the collisional transition rate coefficients, q^X(T, N_e). Rates for l-changing collisions at a fixed principal quantum number n = 2 are much higher than for n-changing collisions, and $n=2\to 3$ collisions typically require N_e ≳ 10⁶ cm⁻³ (e.g., Hummer & Storey 1987). Therefore, n-changing collisions are ignored.

Pengelly & Seaton (1964) have developed a theory to calculate these collisional transition rates. Vrinceanu et al. (2012) have argued that these rates are significantly overestimated. This, however, has been disputed by Storey & Sochi (2015) and Guzmán et al. (2016). Consequently, NEBULAR uses the prescription of Pengelly & Seaton (1964). Accordingly, the hydrogen two-photon emissivity becomes significantly reduced for N_X ≳ 1000 cm⁻³.

2.3.2. The He i Two-photon Continuum

The 1s2s 1S and 1s2s 3S states in He i may also decay by two-photon emission to the ground state (Drake et al. 1969). The decay rates for 1s2s 3S are ignored, as they are about 10 orders of magnitude lower than for the 1S state. The total transition probability is

$$\begin{eqnarray}&&{A}_{2q}=53.1\,{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 10 }$$

for 1s2s 1S. In analogy to Nussbaumer & Schmutz (1984), the spectral dependence tabulated in Drake et al. (1969) is approximated in NEBULAR,

$$\begin{eqnarray}&&A(y)=C[1+\beta {(y-1/2)}^{2}+\gamma {(y-1/2)}^{4}],\end{eqnarray} \tag{ 11 }$$

with y = ν/ν₁₂, C = 145.29 s⁻¹, β = −2.259 and γ = −8.77. This fit is accurate to 0.9% or better over 0.025 ≤ y ≤ 0.975.

The rest of the treatment is the same as outlined in Section 2.3.1, with the notable distinction that no extensive data are available for α^eff(T, N_e). While CLOUDY calculates these parameters internally, NEBULAR relies on the single temperature dependence between 5000–20,000 K listed in Golovatyj et al. (1997). Fortunately, this is not a big drawback: for typical abundances and depending on T, the He i two-photon continuum is 2–6 orders of magnitude below the H i two-photon continuum. Peak emission in f_λ occurs at around 770 nm.

2.4.1. H i Emission Lines

2.4.2. He ii Emission Lines

2.4.3. He i Emission Lines

2.4.4. Other Remarks

The ionization fractions of H⁺, He⁺ and He⁺⁺ (and their number densities, N_X), are calculated for a mixed H/He plasma in collisional ionization equilibrium, solving the Saha equation similar to Rouse (1962). Therein, the algorithm starts with a fixed matter density ρ, and then iteratively obtains N_e; the ionization fractions follow. In the case of NEBULAR, the iteration is skipped because N_e is a fixed parameter. Figure 1 shows these ionization fractions, calculated for a low and a high density mix. Alternatively, user-defined estimates for the ionization fractions can be given to NEBULAR as command-line arguments.

Zoom In Zoom Out Reset image size

A series of consistency checks have been carried out to minimize the number of bugs in the code.

First, data tables of Pengelly (1964), Drake et al. (1969), Hummer & Storey (1987), Storey & Hummer (1995), Ercolano & Storey (2006), Porter et al. (2012), van Hoof et al. (2014) and Storey & Sochi (2015) are used in the code. They were restructured to allow for easier and homogeneous integration in the software, and some of them had to be digitized manually. This is an error-prone process. All (N_e, T) surfaces of the derived data products were inspected in Mathematica's interactive 3D visualization module; typographical errors, and errors during restructuring of the multi-dimensional data tables are immediately recognizable. For example, a typographical mistake was found in Table VI of Pengelly (1964) for the He ii Pickering series, where it should read 0.39 instead of 3.39 (see Appendix B).

Second, the GSL 2D interpolation performed on the (N_e, T) grid has been verified by comparing the numeric results with those given in the literature (e.g., Ercolano & Storey 2006). For the emission line spectra, Storey & Hummer (1995) and Storey & Sochi (2015) have provided dedicated data servers operating on the original data tables; NEBULAR's GSL 2D interpolation reproduces numerically identical values based on the restructured tables.

Third, some figures in the reference publications compare several quantities with each other. For example, Figure 1 in Ercolano & Storey (2006) displays the free-bound continuum emission coefficient for H i, He i and He ii. A similar plot is found in Figure 4.1 of Osterbrock & Ferland (2006) for the combined free-bound and free–free continuum, also showing the two-photon spectrum. It was verified that the corresponding sub-routines in NEBULAR reproduce these figures (see Appendix A.1 for an example).

Lastly, comparison calculations have been performed with CLOUDY for simple gas configurations. While CLOUDY and NEBULAR rely to a large extent on the same atomic data bases, differences do exist. Most noteworthy, for example, are different amplitudes for the helium free-bound continuum, because CLOUDY is a photo-ionization code, whereas NEBULAR simply solves the Saha equation for a collision ionized gas. To compensate for this, the ionization fractions calculated by NEBULAR may be overridden with command line arguments.

NEBULAR is invoked on the command-line and controlled by mandatory and optional arguments (Appendix C., Table 4). If no arguments are given, NEBULAR reports its usage.

There are three mandatory arguments: a wavelength/frequency range including step size, the electron density N_e, and the electron temperature T. The range is assumed to be a wavelength range in Å if the numeric minimum and maximum values are less than 10⁷; otherwise, frequencies [Hz] are assumed. Alternatively, a file containing an ASCII table may be provided, containing the wavelengths for e.g., each pixel of an observed spectrum.

Three example runs with NEBULAR can be found in Appendix A.

NEBULAR produces an ASCII table with the various components and the total spectrum evaluated at each wavelength (or frequency). The results can be readily displayed in TOPCAT (Taylor 2005).

If a wavelength range is given, then the output spectrum will be in wavelength units (j_λ, erg s⁻¹ cm⁻³ Å⁻¹). If frequencies are given, then the spectrum will be returned in frequency units (j_ν, erg s⁻¹ cm⁻³ Hz⁻¹). Internally, all calculations are performed in frequency units. An excerpt of the output is shown in Appendix C., Table 5. Optionally, instead of j = 1/(4π)N_X N_e γ, the atomic emission coefficients γ can be returned (e.g., to access the tabulated data by Ercolano & Storey 2006; Storey & Sochi 2015). For examples, see Appendix A.

NEBULAR also prints various parameters to the header comment section of the output file, such as the derived ionization fractions, the ionic number densities, and the total physical gas density.

Pure hydrogen (helium) spectra can be produced by setting the abundance to 0 (1).

The synthesis of a complete spectrum from 300–30,000 Å with 0.1 Å resolution takes a few seconds on a 2 GHz CPU. This makes NEBULAR useful to compute e.g., the equivalent width of a particular emission line over the nebular continuum (see Figure 2), and to study the relative contributions of the various continuum processes (Figure 4 in Appendix A.2).

Zoom In Zoom Out Reset image size

Usually, the nebular continuum is weak and often undetected in extragalactic observations. Yet it can become significant at far-UV wavelengths because of the strong two-photon emission. If the continuum has been fixed based on optical observations of e.g., the Hβ line, then its contribution to the GALEX NUV and FUV bandpasses may be estimated (the reason this software was developed, Schirmer et al. 2016). Another application would be the observational properties of the Balmer break at 3646 Å, as a function of spectrograph resolution and dispersion (see Appendix A.3, and also Lintott et al. 2009).

The data for Figure 3 can be reproduced with the following two commands:

• 1.  
  >./nebular -r 1e13 6e15 1e11 -t 10000 -n 1 -u 2 -k ‐o out.dat
• 2.  
  > awk '($0!∼/#/) {print $2, $4+$10, $5+$11, $6+$12, $7}' out.dat > plot.dat

Zoom In Zoom Out Reset image size

The first command calculates the spectrum. Option -r specifies the frequency range and step size, needed to get the result in frequency units. Option -u 2 specifies that the computations yield the emission coefficients (γ_fb etc.), i.e., the atomic data tabulated in e.g., Ercolano & Storey (2006), multiplied by frequency ν. Option -k suppresses the calculation of emission lines. A summary of the syntax is given in Appendix C., Table 4.

The second command line is to extract the wavelengths, the sums of free-bound and free–free emission for the three ionic species, and the hydrogen two-photon spectrum, respectively.

Figure 4 shows a model spectrum. The difference to Figure 3 is that NEBULAR multiplied the various emission coefficients by the electronic and ionic densities. A wavelength range (as compared to a frequency range) is chosen as input, and thus the result is in wavelength units (as opposed to frequency units in the previous example).

Zoom In Zoom Out Reset image size

The data for Figure 4 were calculated with

• 1.  
  >./nebular -r 300 12000 1 -t 14000 -n 100

If one were to compare this spectrum with a real observation, then it just has to be multiplied with a constant factor (effectively, integrating over the nebula's volume and folding in the luminosity distance). This would yield the observable spectral flux density, f_λ(λ).

Zoom In Zoom Out Reset image size

If the observed spectrum is on an arbitrary wavelength scale, then NEBULAR can resample the spectrum on the same grid:

• 1.  
  >./nebular -i $\langle$ ASCII wavelength table $\rangle$ -t 14000 -n 100

The data for this plot can be generated by smoothing the model spectrum with a Gaussian kernel, to mimic the finite spectrograph resolution, R = λ/Δλ. The FWHM of the kernel is set to Δλ. The spectral dispersion is controlled by the wavelength step size (third numeric argument for option -r). For the R = 10,000 and R = 2000 examples, dispersions of 0.12 and 0.60 Å pix⁻¹ are used:

• 1.  
  >./nebular -r 3640 3670 0.60 -t 10000 -n 100 -w 1.80 -b -o R2000.dat
• 2.  
  > ./nebular -r 3640 3670 0.12 -t 10000 -n 100 -w 0.36 -b -o R10000.dat

Options -w and -b specify the kernel width and suppress the output of header lines, respectively.