Abstract
We create a photoionization model embedded in the turbulent interstellar medium (ISM) by using the state-of-the-art Messenger Monte Carlo MAPPINGS V code (M3) in conjunction with the CMFGEN stellar atmosphere model. We show that the turbulent ISM causes the inhomogeneity of electron temperature and density within the nebula. The fluctuation in the turbulent ISM creates complex ionization structures seen in nearby nebulae. The inhomogeneous density distribution within the nebula creates a significant scatter on the spatially resolved standard optical diagnostic diagrams, which cannot be represented by the spherical constant-density photoionization model. We analyze the dependence of different optical emission lines on the complexity of nebular geometry, finding that the emission lines residing on the nebular boundary are highly sensitive to the complexity of nebular geometry, while the emission lines produced throughout the nebula are sensitive to the density distribution of the ISM within the nebula. Our fractal photoionization model demonstrates that a complex nebular geometry is required for the accurate modeling of H ii regions and emission-line galaxies, especially for the high-redshift galaxies, where the ISM is highly turbulent based on increasing observational evidence.
Export citation and abstract BibTeX RIS
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
Robust photoionization models are fundamental to understanding galaxy evolution and chemical evolution. Theoretical emission-line models are required to measure the chemical abundance (McGaugh 1991; Kewley & Dopita 2002; Dopita et al. 2006; Kewley et al. 2019), the electron density (Osterbrock 1989; Kewley et al. 2019), and the pressure of the interstellar medium (ISM) (Kewley et al. 2019). With these diagnostic tools, we can analyze the chemical evolution of galaxies (Kobulnicky & Kewley 2004; Tremonti et al. 2004; Zahid et al. 2011; Yuan et al. 2013; Strom et al. 2022), the star formation history of galaxies (Hopkins & Beacom 2006; Madau & Dickinson 2014; Tacchella et al. 2022), and the power sources of galaxies across cosmic time (Kewley et al. 2013; Rigby et al. 2021).
One of the most powerful applications of theoretical photoionization models has been the determination of power sources of galaxies using the strong emission lines available in the optical spectra of galaxies. Baldwin et al. (1981) first proposed the usage of the [N ii]/Hα and [O iii]/Hβ ratios to classify the excitation sources in galaxies, which is known as the Baldwin–Phillips–Terlevich (BPT) diagram. The [N ii]/Hα and [O iii]/Hβ ratios excited by active galactic nuclei (AGNs), shocks, and planetary nebulae are larger than the ratios excited by star formation due to the difference in the hardness of ionizing radiation fields. Veilleux & Osterbrock (1987) derived a semiempirical classification scheme and extend the standard optical diagnostic line ratios to include the [S ii]/Hα and [O i]/Hα ratios.
Kewley et al. (2001) derived the first theoretical classification scheme for local galaxies by computing photoionization and shock models with the MAPPINGS III photoionization code (Binette et al. 1985; Sutherland & Dopita 1993) in conjunction with the stellar population synthesis models (Leitherer et al. 1999). Kewley et al. (2013, 2013) further model the changes of the diagnostic line ratios as a function of cosmic time, finding the [N ii]/Hα and [O iii]/Hβ ratios are likely to be large in high-redshift galaxies due to the extreme ISM conditions at high redshift (Shapley et al. 2005; Steidel et al. 2014).
All previous models for the spectra in star-forming galaxies assume a spherical or plane-parallel geometry. However, spatially resolved spectra of local H ii regions show the complex structures within the nebula. Kinematic and morphology studies reveal the filaments (Gendron-Marsolais et al. 2018; Zavala et al. 2022) and turbulent structures (Rubin et al. 2011) within the nebula
These complex structures are not spherical or plane parallel (Y. F. Jin et al. 2022, in preparation). Snijders et al. (2007) show that a clumpy photoionization model is required to reproduce the extremely high ionized densities and ionization parameters as found in the local starburst galaxies M82 and the Antennae.
In M82 and the Antennae, the young massive H ii regions are studied as candidates for the star-forming regions in high-redshift galaxies through their dense and compact morphology (Smith et al. 2006) and the large ionization parameters (Snijders et al. 2007; Rigby et al. 2011). In particular, the large ionization parameters seen in local extreme star-forming regions appear to be widespread in the current set of high-redshift galaxies observed using near-infrared spectroscopy (Shirazi et al. 2014; Sanders et al. 2016; Papovich et al. 2022).
In this Letter, we present the first geometrical photoionization model using over 30 chemical elements embedded in the turbulent ISM by using M3 (Messenger Monte Carlo MAPPINGS V; Jin et al. 2022), with detailed stellar atmosphere models as the input ionizing source. We present the complex structures of the ionizing radiation field, the electron temperature, the electron density, and the ionization states of different elements within the model. We investigate how the nebular geometry affects the standard optical diagnostic emission lines in the scope of both spatially resolved and integrated spectroscopy. Our results have important implications for studies of emission-line ratios in high-redshift galaxies.
2. The Method
M3 (Jin et al. 2022) is a descendant of the MAPPINGS V photoionization code (Sutherland & Dopita 1993). M3 combines the Monte Carlo radiative transfer technique (Lucy 1999) with MAPPINGS V (Binette et al. 1985), which comprehensively models the microphysics of over 30 chemical elements in the ISM. M3 is designed to create photoionization models with arbitrary geometry in three dimensions by considering both stellar and diffuse ionizing radiation fields.
Our model is set up in a fractal ISM density field (hereafter the fractal model) to mimic the hierarchical structures of turbulence in the ISM (Larson 1981). Following Sutherland & Bicknell (2007), we generate the fractal ISM density field by adopting a standard log-normal density distribution with μ = 1.0 and σ2 = 5.0 and a Kolmogorov spatial structure power law. The average density of the ISM density field is nH = 100 cm−3. The static log-normal density distribution is mathematically well constrained and shows the structure and density fluctuations similar to real turbulent molecular clouds. The fractal ISM density field reproduces the clumps and turbulent structures seen in the real nearby H ii region, like Orion (Rubin et al. 2011). We adopt the solar abundance based on the results from Asplund et al. (2009). In this Letter, we are creating a dust-free photoionization model without dust depletion.
We select the ionizing spectrum from the CMFGEN stellar library (Hillier & Miller 1998, 1999). The CMFGEN library accurately describes the UV stellar spectra by using a non-LTE stellar atmosphere. In this work, we choose an O star with luminosity L = 106 L⊙, temperature Teff = 40,000 K, and gravity lg(g) = 4. The selected spectrum well describes the hard ionizing radiation from massive stars (Simón-Díaz & Stasińska 2008).
Our model is set up in a cube with 129 × 129 × 129 spaxels (spatial pixels). The physical size of the entire computational domain is 43 × 43 × 43 pc3. The size is selected based on the size–density relation for the H ii regions in both star-forming galaxies and starburst galaxies (Hunt & Hirashita 2009). A spatial resolution of 0.25 pc per spaxel guarantees that the substructures of the nebulae can be resolved.
We also produce an equivalent spherical nebula model as a reference to illustrate the geometric effects on emission-line distributions, emission-line fluxes, and ratios. The spherical model is created by MAPPINGS V with an assumption of spherical nebular geometry, in conjunction with the stellar ionizing spectrum from the CMFGEN stellar library. Except for the geometry, the input parameters of the spherical model are the same as the input parameters of the fractal model.
3. Results
3.1. The Radiation-bounded Nebula
Figure 1 presents the 2D map of the integrated Hβ emission-line luminosity of the fractal model. The nebula has a fractal boundary consisting of pixels ionized by both stellar and diffuse photons. The nebula shape is determined by the surrounding ISM clumps. The photoionized region extends farthest in low-density regions around the ionizing source. The low density of the ISM in these regions allows photons to easily pass through. In the dense regions, the ionizing photons are absorbed by the dense clumps, leading to a nebular boundary much closer to the central ionizing source. Beyond the dense nebular boundary, there are regions ionized by diffuse photons, which are scattered by dense clumps. These regions, which are ionized by the diffuse photons rather than the stellar photons, are defined as the diffuse ionized regions.
Figure 1. Upper panel: the projected Hβ luminosity map of the 3D nebula with fractal geometry. Lower panel: the projected distribution of log-normally distributed ISM. The white curve outlines the shape of our modeled nebula. The ionizing source is at the center of the domain (indicated by the white star).
Download figure:
Standard image High-resolution imageThe top-left panel in Figure 2 shows the local ionization parameter distribution across the middle slice of the nebula. The local ionization parameter is defined as
where fν>13.6eV is the flux of the ionizing photons above 13.6 eV through a unit area, NH is the local number density of hydrogen, and c is the speed of light. The ionization parameter indicates the capability of the ionizing radiation field to ionize neutral gas. The average dimensionless ionization parameter is log U = −2.8, within the range of the typical ionization parameter −3.2 < log U < −2.7 for local H ii regions (Dopita et al. 2000) and star-forming galaxies (Moustakas et al. 2010). The diffuse ionized regions have an ionization parameter as low as log U = −5 because these regions are ionized by weak diffuse ionizing photons.
Figure 2. Maps of fundamental parameters across the Z = 0 pc slice of the nebula. Gray contours indicate the distribution of the input of ISM density. We present maps of the dimensionless ionization parameter U, the electron temperature Te , the electron density ne and the distribution of fluxes of the Hα, Hβ, [O iii], [O i], [N ii], and [S ii] emission lines. The nebula is assumed to be at the distance of the Large Magellanic Cloud of 48.5 kpc.
Download figure:
Standard image High-resolution imageThe electron temperature (top central panel in Figure 2) ranges from 5000 to 10,000 K with an average inhomogeneity of 〈Te 〉 = 700 K within the main body of the nebula. The electron temperature increases with radius on average but shows a strong azimuthal asymmetry. The positive temperature gradient is caused by the lack of coolant ions when the ionization parameter declines at a large nebular radius. The temperature is 5000 K at the nebular center and increases to 10,000 K at the edge of the nebula. The diffuse ionized regions have a low electron temperature of 2000 K because the weak ionization field can only weakly heat these regions. The electron temperature shows a spatial variation within the nebula. The high-density clumps have higher electron temperature, and the low-density regions have lower electron temperature.
The electron density distribution (shown in Figure 2, top-right panel) presents a similar log-normal distribution to the neutral gas density. The electron density spans four orders of magnitude from 0.1 to 800 cm−3. The average electron density is ne = 30 cm−3. The inhomogeneity of electron density, which is described by the standard deviation, is 〈ne 〉 = 35 cm−3.
3.2. The Ionization Structure
Figure 2 (lower two panels) presents the internal structure of the six diagnostic emission lines, Hα, Hβ, [O iii]λ5007, [O i]λ6300, [N ii]λ6584, and [S ii]λ λ 6717,31. These emission lines show diverse structures based on their critical densities and ionization potentials.
The Balmer Hα and Hβ lines are produced uniformly throughout the main body of the nebula, extending to the diffuse ionized region. The variation of the Hα and Hβ luminosity is smaller than one magnitude across the entire nebula. The Balmer lines in the diffuse ionized region are fainter than the Balmer lines from the denser ionized gas. On average, the Hα and Hβ luminosity in the diffuse ionized region is 0.8 mag fainter than the average luminosity of the dense ionized region.
The [O iii] emission line is also produced throughout the main body of the nebula, which is ionized by stellar photons. In the diffuse ionized region, the [O iii] emission is fainter than the [O iii] emission from the nebular main body by nine magnitudes. The extremely faint diffuse [O iii] emission is caused by the fact that the hardness of the diffuse ionization field is too low to produce O++ ions, which have the ionization potential of 35.1 eV.
The [O i], [N ii], and [S ii] lines are produced at the outer edge of the nebula. As shown in Figure 3, 99.5% of the [O i] emission, 87.2% of the [N ii] emission, and 87.5% of the [S ii] emission are from the outer 10% radius of the nebula. In the turbulent H ii region model, the [O i], [N ii], and [S ii] lines highlight the "coastline" of the nebula. The fractal nebular boundary extends the perimeter of the outer edge of the nebula, increasing the total luminosity of [O i], [N ii], and [S ii] emission lines compared with a standard spherical model.
Figure 3. Spherical MAPPINGS V model. The spherical model has the same inputs as the fractal model except for the geometry. Here, we present the radial distribution of the [O iii], [O i], [N ii], and [S ii] lines.
Download figure:
Standard image High-resolution imageThe inhomogeneous density distribution of the ISM leads to the substructures of [O i], [N ii], and [S ii] lines in the fractal model. The emission-line fluxes present a clumpy distribution within the nebula, where bright clumps trace high-density regions in ISM. We also find an ionization bar at the center of the nebula crossing from the upper left to the lower right of the z = 0 slice. The substructures seen in the fractal model are similar to the ionization structures seen in the Orion Nebula, where the complex ionization structures coexist with the complex density structures (Rubin et al. 2011).
3.3. The Spatially Resolved Standard Optical Diagnostic Diagrams
We investigate the [N ii]/Hα, [S ii]/Hα, [O iii]/Hβ, [O i]/Hα, and [O iii]/Hβ diagnostics predicted from the fractal model in comparison with the spherical model used in previous theoretical classification schemes for the optical diagnostic diagrams. We collapse the 3D distributions of emission lines into 2D maps to imitate spatially resolved data.
Figure 4 shows the number density distribution of the spaxels from our model on the optical diagnostic diagrams. On all diagnostic diagrams, the spaxels lie on a sequence in the locus of H ii regions bounded by the demarcation between star formation and the harder ionizing sources given in Kewley et al. (2001). The spaxel sequence shows a significant scatter on the diagnostic diagrams with a scatter of 0.8 dex in the log([O iii]/Hβ) ratio, 0.6 dex in the log([N ii]/Hα) ratio, 0.8 dex in the log([S ii]/Hα) ratio, and 1.2 dex in the log([O i]/Hα) ratio. The scatter is defined as the region that includes 16%–84% of the spaxels in the emission-line ratio distributions. The scatter in the diagnostic line ratios is caused by the complex internal structures of the ionization parameter and the hardness of the ionizing radiation field, which vary from spaxel to spaxel, as shown in Figure 2.
Figure 4. The comparison of diagnostic emission-line ratios from the fractal model against the spherical photoionization model. We employ the [N ii]/Hα vs. [O iii]/Hβ (left column), [S ii]/Hα vs. [O iii]/Hβ (middle column), and [O i]/Hα vs. [O iii]/Hβ (right column) diagrams. The background grayscale indicates the number density of spaxels from the collapsed 2D map of the fractal model. The red stars are the spaxels from the imitated long-slit data of the spherical model. The solid line in each panel is the demarcation line given by Kewley et al. (2001) to separate star-forming galaxies/H ii regions from AGNs. The dashed line in the [N ii]/Hα vs. [O iii]/Hβ panel is the demarcation line given by Kauffmann et al. (2003) to separate star-forming galaxies/H ii regions from AGNs.
Download figure:
Standard image High-resolution imageWe collapse the spherical model along the x-axis and the y-axis to produce 1D data to imitate a long-slit observation. Figure 4 shows that the spatially resolved spherical model lies on the spaxel sequence of the fractal model and agrees with the location on the standard optical diagnostic diagrams where the spatially resolved fractal model spaxels are concentrated. The spherical model pixels agree with 19% of the fractal model spaxels on the [N ii]/Hα versus [O iii]/Hβ diagram, 24% of the fractal model spaxels on the [S ii]/Hα versus [O iii]/Hβ diagram, and 20% of the fractal model spaxels on the [O ii]/Hα versus [O iii]/Hβ diagram. However, the spatially resolved fractal model covers a larger range of [N ii]/Hα, [S ii]/Hα, [O i]/Hα, and [O iii]/Hβ ratios than the spatially resolved spherical model.
3.4. Geometric Effects on Total Emission-line Fluxes over a Nebula
Figure 5 gives the deviation between these two models between the spherical and fractal models for the total (integrated) emission-line fluxes for each line.
Figure 5. The comparison of integrated emission-line fluxes and ratios between the fractal model and the spherical model. The emission lines are classified into three categories, which are the Balmer lines (gray shadow), the Volume species (blue shadow), and the Boundary species (red shadow). The black horizontal dashed line indicates no changes in the line fluxes and ratios between the fractal model and the spherical model.
Download figure:
Standard image High-resolution image
The complex geometry makes only a minor change in the fluxes of the Balmer lines of 4% in the Hα flux and 3% in the Hβ flux. This result is expected: Kennicutt & Evans (2012) previously found that the fluxes of Case B recombination Balmer lines are solely correlated with the ionizing luminosity of the central source with limited dependence on the nebular geometry.
The nebular geometry has a pronounced influence on the forbidden lines. We classify the forbidden lines into "Volume species" and "Boundary species" based on their locations within the nebula (Figure 3).
"Volume species" lines are the forbidden lines that are produced throughout the entire nebula and can trace the properties of the overall volume of the ISM within the nebula, [O iii] is an example of a volume-sensitive emission line. In this work, the total flux of the [O iii] line of the fractal nebula model is 31% lower than the flux of the spherical model. The [O iii]/Hβ ratio of the fractal model is also 33% lower than the ratio of the spherical model. The reduced [O iii] flux is caused by the log-normal density distribution where 93% of the nebular volume is occupied by the spaxels with a density lower than the average density of 100 cm−3.
"Boundary species" lines are the ionization lines residing on the boundary of the nebula and are sensitive to the perimeter of the nebular boundary. Complex nebular geometry elongates the nebular boundary, increasing the total flux of "Boundary species" lines. In this work, the [N ii], [S ii], and [O i] lines are the "Boundary species" lines. Compared to the spherical model, the complex geometry of the fractal model increases the total emission-line flux by 29% for the [N ii] line, 161% for the [S ii] line, and 253% for the [O i] line. The total flux of the "Boundary species" is positively correlated with the amount of emission lines concentrated on the edge of the nebula. As shown in Section 3.2, 99.5% of the [O i] emission, 87.5% of the [S ii] emission, and 87.2% of the [N ii] emission is concentrated on the boundary. Therefore, the [O i] line has the largest change in total flux from the spherical to the fractal model among the three lines while the [N ii] has the least change in total flux. Correspondingly, the line ratio is 24% larger for the [N ii]/Hα ratio, 151% larger for the [S ii]/Hα ratio, and 239% larger for the [O i]/Hα ratio for the fractals model, relative to the 1D model.
4. Conclusions
We compute the first fractal-geometry photoionization model for a H ii region with 30 chemical elements to estimate the electron temperature, density, and emission lines throughout the nebula in 3D. Our fractal model shows that the density fluctuation of the turbulent ISM can cause the inhomogeneity of the electron temperature and electron density, and lead to the complex ionization structures seen in the nebula. The fractal geometry more strongly impacts emission-line fluxes of boundary elements [O i], [S ii], and [O iii]. The emission lines of [O iii] and the Balmer lines are only marginally affected by nebular geometry.
We find that the spatially resolved standard optical diagnostic diagrams of the fractal model cannot be reconstructed by an equivalent spherical model, but the average emission-line ratios of these two models have agreement within 0.64 dex for log([O iii]/Hβ) ratio, 0.16 dex for log([N ii]/Hα) ratio, 0.13 dex for log([S ii]/Hα) ratio, and 0.23 dex for log([O i]/Hα) ratio. In conjunction with the complex nebular geometry, the fluctuation of the ISM density within the fractal nebula causes a large scatter in the optical diagnostic diagrams, corresponding to different ionization parameters, densities, and temperature regions across a H ii region. The photoionization model with a uniform density distribution and spherical geometry is insufficient to interpret the spatially resolved observations of the nearby nebula.
Geometric photoionization models are the future in interpreting emission lines from complex ISM conditions, especially for high-redshift galaxies that have evidently turbulent ISM. The ISM turbulence produces the complex H ii region geometry and the density inhomogeneity within nebulae. The perturbations from gas accretion and minor mergers in high-redshift galaxies drive the turbulence of the ISM. Therefore, a more complex nebular geometry driven by turbulence may play a role in the larger [N ii]/Hα ratios seen in high-redshift galaxies.
This research was conducted on the traditional lands of the Ngunnawal and Ngambri people. This research was supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. Y.J. and L.K. gratefully acknowledge the support of L.K.'s ARC Laureate Fellowship (FL150100113).