Messenger Monte-Carlo MAPPINGS V (M^3) -- A self-consistent three-dimensional photoionization code

The Messenger Interface Monte-Carlo Mappings V (M^3) is a photoionization code adopting the fully self-consistent Monte-Carlo radiative transfer technique, which presents a major advance over previous photoionization models with simple geometries. M^3 is designed for modeling nebulae in arbitrary three-dimensional geometries. In this paper, we describe the Monte-Carlo radiative transfer technique and the microphysics implemented in M^3, including the photoionization, collisional ionization, the free-free and free-bound recombination, and two-photon radiation. We put M^3 through the Lexington/Meudon benchmarks to test the reliability of the new code. We apply M^3 to three HII region models with fiducial geometries, demonstrating that M^3 is capable of dealing with nebulae with complex geometries. M^3 is a promising tool for understanding emission-line behavior in the era of SDSS-V/LVM and JWST, which will provide high-quality data of spatially-resolved nearby HII regions and highly turbulent local and high-redshift HII regions.


INTRODUCTION
Modeling nebular emission-line regions is vital for the interpretation of spectroscopic data. Comparing the modeled emission-line spectra with observations can determine the central power source in galaxies, such as star formation, AGN or shocks (Baldwin et al. 1981;Veilleux & Osterbrock 1987;Osterbrock et al. 1992;Kewley et al. 2001;D'Agostino et al. 2019;Byler et al. 2020). Accurate models of emission-line ratios provide the fundamental properties of the interstellar medium (ISM), including the metallicity (Tremonti et al. 2004;Yuan et al. 2013;Yabe et al. 2015;Sanders et al. 2020), the pressure and the electron density of the ionized gas (Kaasinen et al. 2017;Kewley et al. 2019;Harshan et al. 2020;Davies et al. 2021).
Photoionization codes are fundamental tools to model emission-line regions (Ferland et al. 1998;Netzer 1993;Sutherland & Dopita 1993;Ercolano et al. 2003;Moris-set 2006). These models combine atomic data, radiative transfer processes and the physics of the ISM. Ferland (1995) summarize the features of photoionization models. The architectures are similar in most photoionization codes. Nebular models treat the ISM as a plane parallel slab or sphere, dividing the ISM into a series of zones which satisfy ionization and thermal equilibrium (Dopita et al. 2000;Groves et al. 2004). The solution of the radiative transfer equation in each zone is based on the local ionization conditions. The ionization states in each zone are determined by balancing ionization and recombination, including the processes of photoionization, collisional ionization, radiative and dielectronic recombinations. The thermal structure of the nebula is determined by the balance between cooling and heating processes (Sutherland & Dopita 1993). The coefficients for ionization and recombination rates are generated from the fundamental atomic data (See Stasińska 2002, for a review).
The treatment of the diffuse radiation is sophisticated in photoionization models. Approximations are always adopted in photoionization codes because of the com-plexity of calculating the diffuse radiation in radiative transfer processes (Wood et al. 2004). The outwardonly approximation is predominantly implemented in photoionization codes (Ferland et al. 1998;Netzer 1993;Binette et al. 1985), which assumes that the locally produced diffuse radiations follow the outward direction of incident photons. The limitation of the outward-only approximation is that too much energy is trapped in the inner part of photoionized regions, creating an unrealistic temperature structure (Blandford et al. 1990). Alternatively, Harrington (1968) and Rubin et al. (1991) adopted a full treatment of diffuse photons where the diffuse radiations are split into an outward stream and a backward stream, whose contributions are iteratively calculated across the photoionized region. Currently, the full treatment of diffuse radiation only suits simple nebular geometries, like a spherical or axisymmetric geometry.
Simplification of the nebular geometry is another problem in most photoionization codes. The geometry is usually assumed to be spherical or plane-parallel for the analysis of the diagnostic emission-lines (Kobulnicky & Kewley 2004;Levesque et al. 2010) and of the power sources in galaxies (Kewley et al. 2001). However, the realistic structure of Hii regions is too complex to be simply described by spherical or plane-parallel geometry. Detailed observations of the Orion Nebula show a concave structure with a bright dense ionized bar in the foreground of a veil of low-dense ionized gas (Zuckerman 1973;O'Dell et al. 2009). In other nearby nebulae, like the "Pillars of Creation", a large amount of filamentary structures exist (Schneider et al. 2016) due to the joint effects from ISM turbulence and stellar radiation (Gritschneder et al. 2010;Tremblin et al. 2013). These nebular structures complicate the density and temperature structures of Hii regions, altering the fluxes of diagnostic emission-lines, like [S ii]λλ6717,31, [O iii]λ5007, [N ii]λ6583 and Hα .
Monte Carlo radiative transfer (MCRT) provides much promise in modeling realistic photoionized regions in three dimensions. One pronounced feature of MCRT is that it can simulate the genuine radiative transfer process in ISM with an arbitrary geometry. The diffuse radiation in MCRT are fully treated such that the diffuse photons are produced based on the local ISM ionization conditions and are emitted isotropically. The MCRT technique has been successfully applied in some threedimensional photoionization simulations (Ercolano et al. 2003). M 3 is a new self-consistent 3D photoionization code, combining the Monte Carlo radiative transfer technique with the well-known mappings v photoionization code.
The purpose of M 3 is to produce reliable diagnostic emission-lines in the nebulae with arbitrary geometry. mappings v (Binette et al. 1985;Sutherland & Dopita 1993;Sutherland et al. 2018) is a large photoionization code including 30 elements and 80,000 cooling and recombination lines generated from CHIANTI v.8.0 atomic database (Del Zanna et al. 2015). The code selfconsistently calculates nebular cooling and heating processes as well as the complex physics of dust grains.
Compared to other three-dimensional photoionization codes (Ercolano et al. 2003;Wood et al. 2004;Vandenbroucke & Wood 2018), M 3 has a more comprehensive consideration of microphysics in the ISM, including cooling and heating. M 3 inherits engines from mappings v to solve ionization and thermal balance, which are well tested in both observational and theoretical models. This paper is structured as follows. In Section 3, we briefly describe the Monte Carlo radiative transfer technique used in our M 3 code. Sections 4 and 5 describe the microphysics we adopt. In Section 7, we demonstrate the reliability of M 3 by comparing the outputs of benchmark cases between M 3 and other photoionization codes. We present two nebular models with complex geometries produced by M 3 in Section 9. In Section 11, we discuss the capabilities of M 3 . A summary is given in Section 12.

THE ARCHITECTURE OF M 3
The M 3 code is built on two major steps: the Monte Carlo radiative transfer (MCRT), and the thermal and ionization balance calculation. The Monte Carlo radiative transfer is implemented to construct a threedimensional ionization field of the photoionized region, based on which M 3 solves the local balance between ionization and recombination, cooling and heating processes. A set of electron temperature, electron density and ionic fraction of chemical elements are updated, feeding into a new round of Monte Carlo radiative transfer, until the convergence of the electron temperature and the ionic fractions of hydrogen. The detailed working flowchart is presented in Figure 1.

MONTE CARLO RADIATIVE TRANSFER
The main principle of Monte Carlo treatments of radiative transfer is to take photons as the calculation quanta to simulate the local physical processes of ionization and recombination. In the M 3 code, we implement the Monte Carlo radiative transfer method built by Lucy (1999), in which photons of the same frequency, ν, are grouped into energy packets, ε(ν), so that  where n is the number of photons in each energy packet. The use of energy packets rather than individual photons as the calculation quanta improves the computational efficiency (Och et al. 1998;Lucy 1999). To increase the accuracy of the ionizing spectral sampling, we adopt a frequency-dependent energy packet given by where L ν is the luminosity of the source at frequency ν, ∆t is the time of the simulation, and N is the number of energy packets used at each frequency. The direction vectors, (n x ,n y ,n z ), of energy packets are created isotropically in Cartesian coordinates by where α and β are random numbers generated between [0,1]. Each energy packet is then launched from the emitting source, travels through the entire simulated domain and contributes to the local radiation field of the cells it travels through. The mean intensity of the radiation field in each cell is the sum of all passing energy packets, which is given by where J(ν) is the mean intensity of radiation field, N is the number of energy packets passing through the cell, V and l are the volume of the cell and the displacement of each packet. The estimated mean intensity is then fed to the local ionization and recombination processes for calculating the temperature and ionization status.

The trajectories of energy packets
The trajectory of each energy packet must be tracked in order to determine the locations of the absorption events followed by re-emissions of each packet during its journey through the nebula. In a "density bounded" nebula, the trajectories of packets terminate at the boundary of the nebula. In a "radiation bounded" nebula, the packets end their journey after undergoing a specific number of the "absorption -re-emission" loops in the nebula.
We adopt the "cell-by-cell" tracking strategy suggested by Lucy (1999), which follows each energy packet cell by cell along its trajectory in concert to check the occurrence of absorption events. There is an alternative method suggested by Harries & Howarth (1997) that is to first calculate the probability of the occurrence of absorption events as a function of the distance to the central ionizing source. The location of the absorption event is then searched for along the radius of the nebula by comparing a random number against the precalculated probability function. However, this method is less efficient than Lucy (1999) method because it requires searching routines in the computations (Ercolano et al. 2003).

Absorption and re-emission
The absorption is determined by comparing the random optical depth, τ p (ν), with the analytical optical depth, τ l (ν), derived from the physical displacement l. The τ l is given by where κ(ν) and ρ are the absorption coefficient and the volume density of the local cell. If τ l is smaller than the random optical depth τ p , then the packet moves along its initial direction to its next stopping-point in the adjacent cell, where a new τ p is assigned. The stopping-point is the location in each cell where the optical depth is updated. The distance between two stopping points is l. If τ l is larger than τ p , the packet moves the distance of τ p /(κρ), where the absorption event and the follow-up re-emission occur. In this scheme, the τ p (ν) is regulated by a logarithmic formula given by, where p is a random number between [0,1]. The regulation of τ p guarantees the statistical fluxes of the ionizing spectrum follow an exponential decline along the distance from the central source.
The energy packet is re-emitted in situ immediately once absorption occurs. The frequency of the re-emitted energy packet is determined by the probability density function (PDF) derived from the local emissivity distribution. The probability density function is given by, where I ν is the emissivity distribution of the diffuse radiation, ν min and ν max bracket the frequency range of the diffuse ionizing spectrum, and p(ν) denotes the probability of an energy packet created at the frequency ν. The local emissivity includes the recombination emissionlines and the continuum recombination.

Photoionization and Collisional Ionization
In M 3 , the photoionization cross section is calculated using the method suggested by Seaton (1958), with the photoionization data given by Raymond (1979), Osterbrock (1989) and Gould & Jung (1991). We also consider the Auger transition by including the modification factor from Weisheit (1974).
The collisional ionization rates are calculated from the electron energy distribution. The Maxwellian distribution is predominantly adopted in Hii region models where the gas is assumed to be completely thermalized (Spitzer 1941;Draine & Kreisch 2018). Recently, the κ electron energy distribution has been found in interplanetary plasmas (Pierrard & Lazar 2010) and used to account for the observed line ratios in Hii regions (Binette et al. 2012;Nicholls et al. 2012Nicholls et al. , 2013.
The electron energy distribution can be chosen to be a Maxwellian distribution or a κ-distribution in M 3 . Under the assumption of a Maxwellian distribution, the collisional ionization rates are calculated based on the method proposed by Arnaud & Rothenflug (1985) and Younger (1981). This method uses a five parameter fit to the collisional cross section and derives an expression for the integral over a Maxwellian velocity distribution in electron energy. For the κ-distribution, we add a theoretical correction factor to the photoionization rate for each atomic shell (Nicholls et al. 2012(Nicholls et al. , 2013.

Recombination
We include radiative recombination and dielectronic recombination for both hydrogenic and non-hydrogenic ions. The recombination rates for hydrogenic ions and non-hydrogenic ions are calculated separately. We use the power-law fits proposed by Aldrovandi & Pequignot (1973) to calculate the recombination rates for all nonhydrogenic ions.
The recombination rates for hydrogenic ions are treated carefully because hydrogen is the dominant species in the ISM and has a critical impact on nebular properties.
We follow the strategy proposed by Sutherland & Dopita (1993), which provides a temperature-dependent calculation of recombination rates. The Seaton (1959) recombination rates are used when log(T /Z 2 ) < 6.0, where T is the temperature and Z is the nuclear charge number. When log(T /Z 2 ) ≥ 6.0, a "quantum correction" factor (Gaunt 1930) is added to the recombination rates in order to avoid the divergence of the Seaton (1959) calculation in this temperature range.

Continuum Radiation and Line Radiation
The continuum radiation in M 3 consists of three parts: free-free, free-bound and two-photon radiation. We adopt the Gaunt factors given by Gronenschild & Mewe (1978) for free-free and two-photon radiation. The Gaunt factors for free-bound radiation are determined by the expression in Mewe et al. (1986).
The line radiation in M 3 includes resonance lines, forbidden lines, inter-system lines and fine-structure lines. The hydrogen and helium resonance lines are given as the linear combination of the Case A and Case B emissivities. The remaining resonance lines are calculated based on the expressions given by Mewe & Gronenschild (1981). In the current version of M 3 , we assume that the resonance lines have the same cross section as the continuum emission. The forbidden-line emissivities are computed by treating the atoms and ions as a five-level system. Because the forbidden lines are always optically thin, we assume that the forbidden-line emissions carry away energies from the ISM without further interaction. Finally, the inter-system lines and fine-structure lines are treated as the radiation from a two-level system using the method suggested in Binette et al. (1982).

HEATING AND COOLING PROCESSES
M 3 inherits the cooling and heating calculations from mappings v, which are described in detail in Sutherland & Dopita (1993). Briefly, the net cooling function is defined as the difference between the energy lost through cooling processes and the energy gained through heating processes. The collisional line, free-free and twophoton radiations are the major cooling sources. The heating mechanisms consist of photoionization heating, collisional ionization and Compton heating. The radiation from recombination processes may contribute to either heating or cooling effects, depending on the electron energy distribution.

DECOMPOSITION OF SIMULATED DOMAINS
High spatial and frequency resolved Monte-Carlo photoionization models require CPUs with large memory. In order to run photoionization models with large grids, we divide the simulated domain into a series of equal-sized blocks, which are distributed to separate CPU-cores on supercomputers. These blocks contain the same number of cells and the adjacent blocks share the same cell wall. Each block contains the data to compute the local thermal and ionization status independently. The frequency and direction vector of energy packets, are transferred between adjacent blocks through the MPI (Message Passing Interface) message passing standard.

TESTING M 3 WITH THE LEXIN GT ON/M EU DON STANDARD MODELS
We apply our M 3 code to the Lexington/Meudon models. The Lexington/Meudon standard models are a series of artificial cases designed by modellers at the workshops in Meudon, France (Péquignot 1986) and Lexington, Kentucky, USA (Ferland 1995), for the purpose of testing the capability of each code in modelling lowtemperature, high-temperature Hii regions and planetary nebula.
We run three Lexington/Meudon benchmarks, HII20, HII40 and PN150. All the three benchmarks are run in a low-resolution mode with 33 3 cells and a high-resolution mode with 55 3 cells. All the models are run on the computer with a 2.9-3.3GHz CPU. The HII20 model represents the physical conditions in low-temperature Hii regions with 20000 K and the HII40 model represents high-temperature Hii regions with 40000 K. The PN150 model is representative of planetary nebulae which have hard ionization fields compared to normal Hii regions. The details of each model are listed in Table 1. Figure 2 presents the radial profiles of electron temperature and the ionic fraction predicted by M 3 . We reproduce the tight temperature and ionic fraction profiles. The He + fraction scatters largely beyond the He + ionizing front in the HII20 model, which is from the residual ionization caused by charge exchange reactions between hydrogen and helium. However, the value of He + fraction is below 1 per cent, which has no impact on the temperature, electron density and emission-line fluxes.
Predicted emission-line fluxes are listed in Table 2 to  Table 4. We also show the results given by other photoionization codes and calculate median fluxes for each emission-line as a reference. These reference photoionization codes are cloudy (Ferland 1995, hereafter GF), mappings v (Sutherland & Dopita 1993, hereafter RS), mocassin3d (Ercolano et al. 2003, hereafter BE), nebula (Rubin et al. 1991, hereafter RR) and P. Harrington's code (hereafter PH). cloudy and mappings v are two popular one-dimensional photoionization codes using the outward-only approximation on the treatment of diffuse radiation field. The RR and PH codes are chosen because they are the only two codes solving the diffuse radiation transfer by iterative calculations. mo-cassin3d is a Monte Carlo photoionization code developed with a different treatment of microphysics. Compared with mocassin3d, M 3 has a more comprehensive consideration of the ionization and recombination processes in the ISM.
In the HII20 and PN150 models, all emission-line fluxes predicted by M 3 are within 2σ deviation to the reference fluxes. In the low-resolution mode HII40 model, [O i]λ6300 + 6363 and [S ii]λ6716 + 6731 are 2 − 3σ of the reference emission-line fluxes. In the high-resolution mode HII40 model, the Hβ fluxes are 2 − 3σ of reference emission-line fluxes.

THE EFFECT OF SPATIAL RESOLUTION ON EMISSION-LINE FLUXES
We compare the emission-line fluxes between the low-resolution and high-resolution M 3 models. Most emission-line predictions are consistent between the two resolution modes with less than 10 per cent difference in flux.
In contrast, the differences in emission-line fluxes are greater than 10 per cent for  Deharveng et al. 2015), which correspond to the blister, bipolar, spherical and fractal Hii region models respectively. The spherical Hii region can be modeled by previous one-dimensional photoionization codes. We apply M 3 to model three fiducial complex geometries of real nebulae, which are blister, bipolar and irregular Hii regions, to display the capability of M 3 in dealing with complex geometries of nebulae.

Blister HII regions
The blister model is a simplified Hii region model, where half of the Hii region is embedded in the dense molecular cloud and half of the Hii region is expelling into the low-dense clumpy ISM (Tenorio-Tagle 1979; Duronea et al. 2012;Panwar et al. 2020). Figure 3 shows the input neutral hydrogen density distribution through the middle panel (x=0) of the blister model. The ISM density distribution consists of three parts: a low density component with n H = 10 cm −3 , an intermediate density cloud with n H = 500 cm −3 and a high density zone with n H = 1000 cm −3 . The value of 10 cm −3 is consistent with the average density of the hydrogen atom in interstellar space (Brinks 1990). The intermediate density cloud is designed to reproduce the clumpiness of the ISM. The high density zone is designed as a "bowl-like" shape, consistent with the scenario that blister Hii regions are caused by the stars formed at the edge of the cloud (Gendelev & Krumholz 2012).
We select a blackbody with temperature of 40000 K as the simple ionizing source. The ionizing source is placed at the center of the simulated domain, with a total luminosity L tot = 3.1 × 10 39 erg s −1 . The inner radius of nebula is R in = 3×10 18 cm. We adopt the solar abundance set given by Asplund et al. (2009) (hereafter AS09) as the ISM chemical abundance. Figure 4 presents the three-dimensional shape of the blister Hii region model. The blister model produces an azimuthally-symmetric nebula consisting of two components. The upper component is the emission from a partially ionized clump and the lower component is the emission from the high density ISM. M 3 successfully reproduces the ionization cone and the diffused ionized gas behind the intermediate gas clump. The high density and low density gaseous components represent dense natal molecular clouds and the ionized bubble created by ionization front, stellar wind and photon pressure (Gendelev & Krumholz 2012). The intermediate density gaseous component represents clumps in the ionized bubble.
We further display the line-of-sight emission-line maps of the blister model in Figure 4 to mimic imaging observations. Similar to spherical nebula models, the Hii region observed along the x-axis has a spherical geometry. However, the integrated emission-line map shows two separate components in views along the y-axis and z-axis.
We  Figure 5 shows the slices of electron temperature, electron density and H-ionizing photon flux across planes of x=0, y=0 and z=0. The electron temperature, the electron density and the H-ionizing photon flux are uniformly distributed within the main body of the nebula. The main body of the nebula has an average electron temperature of 4000 K and an average electron density of 10 cm −3 . The intermediate density cloud has a hotter average electron temperature of 6000 K and a higher average electron density of 450 cm −3 than the remaining part of the nebula. The average electron temperature in the diffuse ionized gas is 1000 K cooler than the main body of the nebula. The average electron density in the diffuse ionized gas is 10 cm −3 similar to the main body of the nebula. Figure 6 shows the slices of ionic fractions of H + , He + and O ++ across planes of x=0, y=0 and z=0. The H + and He + are dominant ionic species within nebulae and the O ++ is the major coolant species. The ionic fractions of H + , He + and O ++ are uniformly distributed across the nebula.

Bipolar HII regions
The bipolar Hii region model is composed of two ionized lobes located perpendicularly to a dense molecular cloud (Samal et al. 2018). Three-dimensional simulations suggest that the bipolar structure is the consequence of a star evolving in a sheet-like molecular cloud (Bodenheimer et al. 1979;Wareing et al. 2017). Figure 7 shows the initial condition of neutral hydrogen density distribution through the middle panel (x=0) of the bipolar Hii region model. The ISM density distribution consists two lobes with n H = 100 cm −3 oriented  perpendicular to a sheet-like high density cloud with n H = 10000 cm −3 . The density of 100 cm −3 within two lobes is selected to match with the typical density of local nebula in star-forming galaxies (Kewley et al. 2001;Levesque et al. 2010). The density of the high density cloud is selected based on data that indicate that the sheet-like clouds have density of 10 4 cm −3 (Deharveng et al. 2015). The central ionizing source is selected as a blackbody with temperature of 40000 K and a total luminosity L tot = 3.1 × 10 39 erg s −1 . The inner radius of nebula is R in = 3 × 10 18 cm. The ISM chemical abundance is the AS09 solar abundance. The bubble structures in the bipolar model exist in most Hii regions (Churchwell et al. 2007;Deharveng et al. 2015). Bipolar structures may appear as ring-like structures if the nebulae are observed in particular viewing angles (Anderson et al. 2011). The 3D visualization of our bipolar model presents ring-like structures similar to real observations. Figure 9 shows the slices of electron temperature, electron density and H-ionizing photon flux across planes of x=0, y=0 and z=0. The radial profiles of electron temperature, electron density and H-ionizing photon flux in the bipolar nebula are similar to the profiles of the spherical model, where the electron temperature, the electron density and the H-ionizing photon flux have flat gradients within the nebula. The average electron density is 108 cm −3 and the average electron temperature is 6000 K. The electron density is reduced and the electron temperature becomes cooler at the edge of the nebula because the ionizing photons are totally absorbed. Figure 10 shows the slices of ionic fractions of H + , He + and O ++ across planes of x=0, y=0 and z=0. The ionic fractions of H + and He + are uniformly distributed in the nebula. In contrast, the ionic fraction of O ++ decreases along with the radius of the nebula.

Fractal Geometry HII regions
The fractal Hii region model represents an Hii region evolving in the turbulent ISM (Medina et al. 2014) and developing self-similar structures (Zuckerman 1973;Rubin et al. 2011;Arthur et al. 2016;O'Dell et al. 2017). Both observations (Wisnioski et al. 2015) and simulations (Pillepich et al. 2019) reveal that the turbulent motion dominates the ISM kinematics especially in highredshift galaxies. Figure 11 shows the neutral hydrogen density distribution in the fractal model. Following the Kolmogorov's theory (Kolmogorov 1941), the power spectrum of the turbulence is in the format of E(k) ∝ k −5/3 , where k is the wavenumber k ∼ 1/r. The column density probability distribution function (PDF) of the ISM is designed as the log-normal distribution because of the hierarchi- cal structures of the turbulent ISM (Larson 1981). The mean value of the column density is 100 cm −3 . We use AS09 for the ISM chemical abundances. The central ionizing source is a blackbody with temperature of 40000 K and a total luminosity L tot = 3.1 × 10 39 erg s −1 . The inner radius of nebula is R in = 3 × 10 18 cm. Figure 12 displays the shape of the fractal Hii region model in three-dimensions. The modeled Hii region has an inhomogeneous geometry which is caused by the fractal density distribution of the ISM.
Theories suggest that a turbulent model best describes the ISM which has a fractal density distribution (Federrath et al. 2009   which neither spherical models nor plane-parallel models can produce. We also present the line-of-sight integrated emissionline maps of the fractal model. The Hii region morphology changes in different line-of-sight views. In each line-of-sight direction, the Hβ and [O iii] share the similar spatial distribution, which are concentrated around the center of the nebula. Compared with the Hβ and [O iii], the [N ii] is less concentrated to the central source because the high local ionization parameter ionizes N + to higher ionization levels. Figure 13 shows the slices of electron temperature, electron density and H-ionizing photon flux across planes of x=0, y=0 and z=0. The electron temperature and the H-ionizing photon flux are uniformly distributed across the nebula. In contrast, the electron density distribution shows a large fluctuation, where the high density ISM has the high electron density and the low density ISM has the low electron density. The average electron density is 150 cm −3 with the standard deviation of 90 cm −3 . Figure 14 shows the slices of ionic fractions of H + , He + and O ++ across planes of x=0, y=0 and z=0. The ionic fractions of H + and He + are uniformly distributed across the model. The ionic fraction of O ++ decreases as a function of the nebular radius.

COMPARISON BETWEEN MODELS WITH COMPLEX GEOMETRIES AND SPHERICAL MODELS
We create a corresponding spherical model for the blister, bipolar and fractal nebula models respectively. Ex-  cept for the geometry, the spherical model has the same initial conditions as the geometric models, including the ionizing source, the ISM abundance, the total mass of nebula and the average ISM density. The corresponding spherical models to the blister and bipolar Hii regions are radiation-bounded because the radiation field is absorbed before reaching the edge of ISM. The corresponding spherical model to the fractal Hii region is density-bounded, because the ISM density inhomogeneity may shorten the radiation field. Figure 15 presents the comparison of emission-line distributions between the geometric models and the spheri-    (Rubin et al. 2011;Peimbert 2019). Detailed observations of θ 1 Ori C in the Orion Nebula show that the electron temperature has an increasing trend with the distance from the ionizing star. In current constant density photoionization mod-els, the electron temperature becomes hotter towards the edge of Hii regions when the metallicity of nebulae 12 + log(O/H) ≥ 8.5 .
The mechanisms leading to the electron temperature fluctuations are still under debate. Multiple mechanisms have been proposed, including the turbulence and shocks in the ISM (Peimbert et al. 1991;O'Dell et al. 2015), and the stellar winds produced by central ionizing sources (Gonzalez-Delgado et al. 1994) Here we propose that the nebular geometry is one, but not the only cause of the electron temperature fluctuation. The degree of electron temperature fluctuation increases with the complexity of nebular geometry. The electron temperature in the blister and bipolar Hii region models has a flat gradient because the hydrogen density is uniformly distributed within ionized bubbles. In the blister Hii region model, the intermediate density clump is hotter than the average electron temperature of the nebula because the clump is denser than the average density of the nebula. The fractal Hii region model has a more significant electron temperature fluctuation than the blister and bipolar Hii region models, because the self-hierarchical structure in the fractal model is the most complex geometry among the blister, bipolar and fractal models.
Nebula models allowing the electron temperature fluctuation offers great promise in interpreting emission-line spectra. The spectra of nearby or distant galaxies captured by fixed-size apertures combine the light from an ensemble of Hii regions with complex electron temperature structures. The electron temperature fluctuation is Figure 16. Comparison of integrated emission-line fluxes (left) and emission-line ratios (right) between the models with complex geometries and the spherical models. The horizon dashed line indicates that the integrated fluxes or emission-line ratios from models with complex geometry and from spherical models are the same. a potential cause of the discrepancy between the metallicity determined by recombination lines and those determined with Auroral lines (Peimbert 2019). The turbulent ISM in high-redshift galaxies increase the electron temperature fluctuations of nebulae. Three-dimensional nebula models with complex electron temperature structures provide more realistic predictions of emission-lines than those models with constant temperature, density or pressure assumptions.

Density Structure of HII regions
Realistic nebulae have significant density fluctuations within local Hii regions (Pérez et al. 2001;Simpson et al. 2004;McLeod et al. 2016;Peimbert 2019). Spatiallyresolved measurements show negative density gradients in some local Hii regions (Kurtz 2002;Phillips 2007;Rubin et al. 2011) and flat density gradients in other compact Hii regions (García-Benito et al. 2010;Ramos-Larios et al. 2010). Detailed observations of the Orion Nebula find that the turbulence within the nebula leads to the complex density structures (Arthur et al. 2016;O'Dell et al. 2017;Kewley et al. 2019).
We present diverse electron density distributions in three fiducial Hii region models. The bipolar Hii region model presents a flat electron density gradient because the density distribution of hydrogen is uniform within the bipolar bubbles. The blister and fractal Hii regions present significant electron density fluctuations, which are caused by the inhomogeneity of the ISM density within the nebulae.

Emission-Line Predictions
The and Hβ. The inhomogeneity of ISM density changes the spatial distributions of emission-lines within nebulae, which cannot be reproduced by nebula models with spherical or plane parallel geometries. We will investigate the detailed impact of nebular geometry on the optical emission-lines in a forthcoming paper (Jin et al. in prep).
Oversimplified nebular geometries in emission-line modelings lead to the debate on how the nebulae are bounded. Spherical and plane parallel nebula models assume that Hii regions are simply radiation-bounded, in order to reproduce the observed emission-line fluxes in nearby galaxies (Kewley et al. 2001). However, Nakajima et al. (2013) proposed that the Hii regions in high-redshift galaxies are likely to be density-bounded, where the gas is insufficient to absorb the entire ionizing photons from the star. Detailed studies of nearby Hii regions suggest that the radiation-bounded and density-bounded conditions coexist within a single nebula (Pellegrini et al. 2012). M 3 can produce emission-lines in Hii region models with complex geometries, allowing a mixture of radiation-bounded and density-bounded cases. Among the three fiducial nebula models, the bipolar Hii region model and the fractal Hii region model are simply radiation-bounded cases. The blister Hii region model is a mixture of radiation-bounded and density-bounded nebula, where the nebula is radiation-bounded in the high density zone and the intermediate density cloud, and the nebula is density-bounded in the low density gaseous component.

CONCLUSION
M 3 is a photoionization code designed for modeling the Hii regions with arbitrary three-dimensional geome-tries. M 3 incorporates the Monte Carlo radiative transfer technique with the complete ISM microphysics implemented in mappings v code, producing realistic threedimensional ionization and thermal structures within nebulae. The accurate cooling functions in M 3 promise reliable predictions of the emission-line fluxes.
We put M 3 successfully through the Lexington/Meudon benchmarks test, which is a series of artificial photoionized models accounting for the physical conditions of various types of Hii regions and planetary nebulae. The emission-line fluxes predicted by M 3 are consistent with the fluxes produced by the reference photoionization codes. We run each Lexington/Meudon benchmark with a high-resolution mode and a low-resolution mode, finding the spatial resolution effect is pronounced for emission-line fluxes produced at the edge of Hii regions.
We create three fiducial Hii region models with complex geometries, which are the blister geometry, the bipolar geometry and the fractal geometry. We find that: • In the blister Hii region model, the high-density clump is partially ionized. The high-density clumpy structure has hotter electron temperature and higher electron density than the low-density gas in the nebula. The diffuse ionized gas is partially ionized and cooler than the average electron temperature of the nebula.
• The bipolar Hii region model has the similar radial profiles to the spherical model, in terms of the electron temperature, the electron density, the ionizing photon flux, the ionic fraction and the emission-line intensities. However, the bipolar Hii region has smaller volume size than the spherical model, reducing the integrated emission-line fluxes.
• The fractal model has the most complex geometry among the three fiducial Hii region models. Neither the internal ionization and thermal structures nor the integrated emission-line fluxes can be reproduced by simple spherical photoioinzation models. The inhomogeneity of ISM density causes the fluctuation of the electron temperature and the electron density. The twisted nebular boundary of the fractal model increase the boundary emissionline species.
We demonstrate that M 3 is a promising tool for interpreting nebular emission-line behaviors in the era of JWST and the upcoming local integral-field unit surveys, like SDSS-V/LVM (the Local Volume Mapper).
We are grateful for the valuable comments on this work by an anonymous referee that improved the scientific outcome and quality of the paper. This research was conducted on the traditional lands of the Ngunnawal and Ngambri people. YFJ is grateful to Brent Groves, Christoph Federrath, Emily Wisnioski and David Nicholls for useful discussions. 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. L.J.K. gratefully acknowledges the support of an ARC Laureate Fellowship (FL150100113).