NEBULAR: A simple synthesis code for the hydrogen and helium nebular spectrum

NEBULAR is a lightweight code to synthesize the spectrum of an ideal, mixed hydrogen and helium gas in ionization equilibrium, over a useful range of densities, temperatures and wavelengths. Free-free, free-bound and two-photon continua are included as well as parts of the HI, HeI and HeII line series. NEBULAR interpolates over publicly available data tables; it can be used to easily extract information from these tables without prior knowledge about their data structure. The resulting spectra can be used to e.g. determine equivalent line widths, constrain the contribution of the nebular continuum to a bandpass, and for educational purposes. NEBULAR can resample the spectrum on a user-defined wavelength grid for direct comparison with an observed spectrum; however, it can not be used to fit an observed spectrum.


Introduction
Astrophysical nebulae have complex spectra. Photoionization codes, such as CLOUDY , calculate realistic spectra of nebulae for a wide range of physical parameters of the gas and the ionizing sources near or within them. These codes may have complex setups and significant execution times because of their ambitious scope.
NEBULAR is not a contender in this area, because it does not fit an observed spectrum, nor does it return physical parameters. However, it synthesizes an accurate and fairly complete spectrum of an ideal, mixed hydrogen and helium gas in ionization equilibrium, over a wide range of temperatures, densities and wavelengths. Hydrogen and helium are the most abundant elements in the Universe, and thus form a central part in the spectral analysis of gaseous nebulae. NEBULAR provides easy and quick access to the most recent atomic computations of these elements. Because it is lightweight, it could be used in larger analysis packages. It is expedient for physical purposes and astronomy alike, whenever the general properties of a spec-trum are investigated. For example, it can take the arbitrary wavelength grid of an observed real spectrum as input and compute a model spectrum on the same grid for direct comparison. Because of the simple installation and usage, it can also be used for teaching purposes.
NEBULAR is open source, licensed under the GNU Public License (GPL) and available on Github 1 ; this paper refers to version 1.2. NEBULAR is written in C++ and links against common standard libraries, facilitating cross-platform usability and installation. This paper is organized as follows. Methods and relevant physical relations are recapped in Section 2, together with explanations of the code's limitations. Section 3 gives an overview of software verification, usage, the output format, and a few applications. A summary is given in Section 4. Three example runs are presented in the Appendix, together with further ancillary tables. All wavelengths in this paper are vacuum wavelengths, and all logarithms refer to base 10.

Methods and limitations of the code
The following is a summary of the physical processes involved.
The combined hydrogen and helium spectrum is synthesized from four different components: the free-bound, free-free and two-photon continua, and part of the H i, He i and He ii emission line series. The respective rescaled emission coefficients j are written as Here, N X and N e are the ion and electron number densities in cm −3 , respectively, and γ is the emission coefficient (in erg s −1 cm 3 Hz −1 ; apart from j nn , which has no frequency dependence). The sources of the input data, as well as the limitations of the code, are reviewed below.

Free-bound continuum
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 2 ≤ T ≤ 10 5 K with a step size of ∆ log 10 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 continuum
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 coef-ficients are Brown & Mathews 1970). Here, e, m e , h, and k B are the electron charge, electron mass, Planck constant and Boltzmann constant, respectively. hν 0 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 −10 statC, with 1 statC = 1 cm 3/2 g 1/2 s −1 . 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 g ff (γ 2 , u) for Maxwell-Boltzmann electron distributions, where γ 2 = Z 2 hν 0 /(k B T ) and u = hν/(k B T ).
The data for g ff (γ 2 , u) are available on a loglog 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).

H i and He ii
The 2 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 2 S level by emission of a single photon 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 2 S → 1 2 S two-photon decay is (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 (Hi- 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 (Sect. 2.1) below 2000Å .
The three factors are discussed in the following.
Effective recombination coefficient α eff is the effective recombination coefficient to the 2 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 the Appendix (Tables 4 and 5). For convenience, these tables contain the coefficients for the 2 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 α eff 2 2 S (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 α eff 2 2 S (T, N e ) at these densities and temperatures is weak.

Frequency dependence g(ν)
The frequency dependence of the two-photon spectrum is absorbed in Here, ν 12 is the Lyα frequency of the respective hydrogenic ion. Nussbaumer & Schmutz (1984) have presented an analytic fit for the frequency dependent transition probability, A(ν/ν 12 ). It is accurate to within 1% over nearly the full spectral range, Collisional de-population of the 2 2 S level The 2 2 S level can be de-populated by means of angular momentum changing (l-changing) collisions, 2 2 S → 2 2 P, inflicted by other ions and electrons. The 2 2 P level may then decay to the 1 2 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 → 3 collisions typically require N e 10 6 cm −3 (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 twophoton emissivity becomes significantly reduced for N X 1000 cm −3 .

The He i two-photon continuum
The 1s2s 1 S and 1s2s 3 S states in He i may also decay by two-photon emission to the ground state (Drake et al. 1969). The decay rates for 1s2s 3 S are ignored, as they are about 10 orders of magnitude lower than for the 1 S state. The total transition probability is for 1s2s 1 S. In analogy to Nussbaumer & Schmutz (1984), the spectral dependence tabulated in Drake et al. (1969) is approximated in NEBULAR, with y = ν/ν 12 , C = 145.29 s −1 , β = −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 Sect. 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 − 20000 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.

H i emission lines
The H i emission line spectrum for Case B (Lyman-thick nebulae) is built from Storey & Sochi (2015) and includes all line transitions from upper levels n ≤ 99 onto lower levels n ≤ 8. The latter was chosen in NEBULAR so that the α line of a j nn series has a wavelength shorter than 36 µm, i.e. it still falls within the free-bound continuum calculated by Ercolano & Storey (2006). The parameter space covers log N e = 2 − 6 and log T = 2.0 − 4.4. The Lyman series (n = 1) is excluded for Case B.
Linear extrapolation is used from log N e = 2 to 0, and has been verified by comparison with Pengelly (1964) for log N e = 0 and T = 1250 − 80000 K. It is accurate within 2.5% for T = 2500− 40000 K. For T = 1250 K and T = 80000 K, the errors are 4% and 7%, respectively. Calculations of the H i line spectrum are limited to upper levels of n ≤ 25 (Storey & Hummer 1995) if extrapolation is used, and are skipped for T > 50000 K.
For Case A (Lyman-thin nebulae), the emission coefficients are taken from Storey & Hummer (1995). The same limitations hold as outlined above; the Lyman series (n = 1) is included .

He ii emission lines
The He ii spectrum for Case A and Case B is based on Storey & Hummer (1995), with a common upper level of n ≤ 25. Treatment is the same as outlined for H i, with some differences: First, transitions onto lower levels n ≤ 14 are included. Second, the temperature range is limited to 500 ≤ T ≤ 10 5 K. Third, calculations are skipped if T > 10 5 K.

He i emission lines
The data from Porter et al. (2013) are used (corrigendum to Porter et al. 2012) to construct the He i spectrum for T = 5000 − 25000 K, log N e = 1 − 14, and Case B. In total, 44 singlet and triplet lines are available between 2946 − 21623Å. Storey & Hummer (1995) and Storey & Sochi (2015) have provided data servers to access and interpolate their data tables. These servers are not directly useful for NEBULAR because of their interactive nature. They were used to extract the data on the pre-defined (log N e , log T ) grid points. As the original grids are not strictly regular, Mathematica was used to build a regular grid, enabling bicubic 2D GSL interpolation for NEBULAR. Within rounding errors, NEBULAR yields the same results as the original data servers.

Other remarks
Recently, it has been discussed that the energy distributions of free electrons may not follow a purely thermal Maxwell-Boltzmann distribution (e.g. Storey & Sochi 2015;Zhang et al. 2016). At this point, NEBULAR adopts Maxwell-Boltzmann distributions throughout.

Ionization fractions
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.

Consistency checks
A series of consistency checks have been carried out to minimize the number of bugs in the code.
First, data tables of Pengelly (1964) 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 multidimensional 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, fig. 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 fig. 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.

Command-line arguments
NEBULAR is invoked on the command-line and controlled by mandatory and optional arguments ( Table 2). 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 7 ; 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.

Output
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 −1 cm −3Å−1 ). If frequencies are given, then the spectrum will be returned in frequency units (j ν , erg s −1 cm −3 Hz −1 ). Internally, all calculations are performed in frequency units. An excerpt of the output is shown in Table 3. 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).

Possible applications
The synthesis of a full spectrum from 300 − 30000Å with 0.1Å resolution takes a few seconds on a 2 GHz CPU. It is straight forward to compute e.g. the equivalent width of a particular emission line over the nebular continuum (Fig. 2), and to study the relative contributions of the various continuum processes (Fig. 4 in Appendix A.2).
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 Fig. 5, and also Lintott et al. 2009).

Summary
NEBULAR is a lightweight code that synthesizes the spectrum of an ideal H/He gas mix in (collision) ionization equilibrium over a wide range of densities, temperatures and wavelengths. It relies on recent atomic computations by others. It's strengths are mainly on the educational side, to study the response of the spectrum and its various constituents to changes in electron temperature and density. Nonetheless, NEBULAR also has scientific applications, e.g. to determine equivalent line widths and to extrapolate fluxes in unobserved wave bands. To infer accurate physical parameters of a real, astrophysical nebula, full photoionization codes such as CLOUDY must be used.
This 2-column preprint was prepared with the AAS L A T E X macros v5.2.

A. NEBULAR example runs
This Appendix contains some example figures and the NEBULAR syntax used to compute the data.
A.1. Extracting atomic data -emission coefficients  Osterbrock & Ferland (2006). Small differences between our data and that of Osterbrock & Ferland (2006) are attributed to the use of recent data from Ercolano & Storey (2006) and van Hoof et al. (2014).
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 Fig. 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). The data for Figure 4 were calculated with > ./nebular -r 300 12000 1 -t 14000 -n 100

A.2. Calculating a model spectrum
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 λ (λ).
If the observed spectrum is on an arbitrary wavelength scale, then NEBULAR can resample the spectrum on the same grid: > ./nebular -i <ASCII wavelength table> -t 14000 -n 100 A.3. Visualizing the pseudo-continuum near the Balmer break Fig. 5.-The hydrogen Balmer break at 3646Å, as seen by two model spectrographs with resolutions R = 2000 and R = 10000, respectively. The confluent high-order Balmer series forms a pseudo-continuum, masking the Balmer break. The broad minimum between 3647Å and 3649Å is artificial, because upper levels higher than n = 99 are not available in NEBULAR.
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). Options -w and -b specify the kernel width and suppress the output of header lines, respectively.

B. Corrections to the literature
This is a list of errors (mostly of typographical nature) found in the literature: • In Table VI of Pengelly (1964) for the He ii Pickering series, for n = 19, t = 1 and Case B, it should read 0.39 instead of 3.39.
• Pengelly & Seaton (1964) erroneously have a factor of π 2 in the denominator for the Debye radius (their eq. 33). Correct would be a single factor π. Subsequent derivations and numeric evaluations, however, are correctly calculated.
• In their eq. (3.40), Brocklehurst (1971) also have a factor of π 2 in the denominator for the Debye radius, and carry it over in their calculations of the collision rates. Spitzer (1956) is cited for the equation, yet in that original work the dependence on π is linear. CLOUDY (v13.03) has the same wrong dependence on π 2 in the calculation of the collision rates, a bug that has been fixed in v15.
• In Table 4.2 of Osterbrock & Ferland (2006) the H i Paschen-line series, the values tabulated for j Pδ /j Hβ actually refer to j P /j Hβ . This table has been derived from Pengelly (1964), where this error is absent.
• In Table 4.3 of Osterbrock & Ferland (2006) the relative intensities with respect to He iiλ4686 are mistakenly a factor of 10 too low for n → 2 and T = 5000, 10000 and 20000 K. This table has been reproduced from Storey & Hummer (1995), where this error is absent. suppresses emission lines if given (faster execution) -u 0, 1, or 2 0: returns j = 1/(4π) N X N e γ (model spectra) 1: returns γ (atomic data) 2: returns νγ Table 3: Resulting output of NEBULAR (abridged). The first two columns contain frequencies and wavelengths. The third column is the total spectrum, either f ν or f λ . The remaining columns give the free-bound, twophoton, free-free and emission-line components for H i, respectively (data for He i and He ii are not shown).