Theoretically Modelling Photoionized Regions with Fractal Geometry in Three Dimensions

We create a photoionization model embedded in the turbulent ISM by using the state-of-the-art Messenger Monte-Carlo MAPPINGS~V code (M$^3$) 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 accurate modeling of HII regions and emission-line galaxies, especially for the high-redshift galaxies, where the ISM is highly turbulent based on the increasing observational evidence.

One of the most powerful applications of theoretical photoionization models has been the determination of power sources of galaxies using strong emission-lines available in the optical spectra of galaxies. Baldwin et al. (1981) first proposed 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 (AGN), 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 semi-empirical 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 MAP-PINGS 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. (2013a,b) 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 in high-redshift (Steidel et al. 2014;Shapley et al. 2005).
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 planeparallel (Jin et al. in prep). 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 of 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 M 3 (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.  (Sutherland & Dopita 1993). M 3 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. M 3 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 (thereafter the f ractal model) to mimic the hierarchical structures of turbulence in ISM (Larson 1981). Follow-ing 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 n H =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 HII 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 CMF-GEN library accurately describes the UV stellar spectra by using non-LTE stellar atmosphere. In this work, we choose an O-star with the luminosity L=10 6 L ⊙ , the temperature T eff =40000 K and the 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 pc 3 . 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 sub-structures 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. 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 furthest 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.

The Radiation Bounded Nebula
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, N H is the local number density of hydrogen and c is the speed of light. The ionization parameter indicates the capability of ionizing radiation field to ionize neutral gas. The average dimensionless ionization parameter is logU =-2.8, within the range of the typical ionization parameter -3.2<logU <-2.7 for local Hii regions (Dopita et al. 2000) and starforming galaxies (Moustakas et al. 2010). The diffuse ionized regions have an ionization parameter as low as logU =-5 because these regions are ionized by weak diffuse ionizing photons. The electron temperature (top central panel in Figure 2) ranges from 5000 K to 10000 K with average inhomogeneity of T e =700 K within the main body of the nebula. The electron temperature increases with radius on average but shows strong azimuthal asymmetry. The positive temperature gradient is caused by the lack of coolant ions when the ionization parameter declines at large nebular radius. The temperature is 5000 K at the nebular center and increases to 10000 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 as the neutral gas density. The electron density spans four orders of magnitude from 0.1 cm −3 to 800 cm −3 . The average electron density is n e =30 cm −3 . The inhomogeneity of electron density, which is described by the standard deviation, is n e =35 cm −3 . 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 magnitude fainter than the average luminosity of the dense ionized region.

The Ionization Structure
The [O iii] emission-line is also produced throughout the main body of 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 diffuse ionization field is too low to produce O ++ ions, which has the ionization potential of 35.1 eV.
The  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 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% to 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 radia-    Kauffmann et al. (2003) to separate star-forming galaxies/H ii regions from AGN. tion field, which vary from spaxel to spaxel, as shown in Figure 2.
We collapse the spherical model along the x-axis and the y-axis to produce a 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 has agreement with the location on the standard optical diagnostic diagrams where the spatiallyresolved fractal model spaxels are concentrated.

Geometric Effects On Total Emission-Lines
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.
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 "V olume species" and "Boundary species" based on their locations within nebula (Figure 3).
"V olume species" lines are the forbidden lines which 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 are occupied by the spaxels with density lower than the average density of 100 cm −3 .
"Boundary species" lines are the ionization lines residing on the boundary of 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. 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 parameter, density and temperature regions across an HII region. The photoionization model with uniform density distribution and spherical geometry is insufficient to interpret the spatially-resolved observations of nearby nebula.
Geometric photoionization models are the future in interpreting emission-lines from complex ISM conditions, especially for high-redshift galaxies which have the 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. YJ and LK gratefully acknowledge the support of Lisa Kewley's ARC Laureate Fellowship (FL150100113).