Optical and JWST Mid-IR Emission Line Diagnostics for Simultaneous IMBH and Stellar Excitation in z~0 Dwarf Galaxies

Current observational facilities have yet to conclusively detect $10^3 - 10^4 M_{\odot}$ intermediate mass black holes (IMBHs) that fill in the evolutionary gap between early universe seed black holes and $z \sim 0$ supermassive black holes. Dwarf galaxies present an opportunity to reveal active IMBHs amidst persistent star formation. We introduce photoionization simulations tailored to address key physical uncertainties: coincident vs. non-coincident mixing of IMBH and starlight excitation, open vs. closed surrounding gas cloud geometries, and different AGN SED shapes. We examine possible AGN emission line diagnostics in the optical and mid-IR, and find that the diagnostics are often degenerate with respect to the investigated physical uncertainties. In spite of these setbacks, and in contrast to recent work, we are able to show that [O III]/H$\beta$ typically remains bright for dwarf AGN powered by IMBHs down to $10^3 M_{\odot}$. Dwarf AGN are predicted to have inconsistent star-forming and Seyfert/LINER classifications using the most common optical diagnostics. In the mid-IR, [O IV] 25.9$\mu$m and [Ar II] 6.98$\mu$m are less sensitive to physical uncertainties than are optical diagnostics. Based on these emission lines, we provide several mid-IR emission line diagnostic diagrams with demarcations for separating starbursts and AGN with varying levels of activity. The diagrams are valid over a wide range of ionization parameters and metallicities out to $z\sim0.1$, so will prove useful for future JWST observations of local dwarf AGN in the search for IMBHs. We make our photoionization simulation suite freely available.


INTRODUCTION
The occupation fraction of supermassive black holes (SMBHs) in massive galaxies is near unity (Magorrian et al. 1998). While LIGO has detected stellar mass black holes (≤ 10 2 M ) resulting from compact object mergers, intermediate black holes (IMBHs) remain elusive in the 10 2 M -10 5 M range . Detecting black holes at the low end of this range would provide a crucial link between early universe black hole seeds and local SMBHs.
Dwarf galaxies show promise in the search for IMBHs based on the M * − M BH relation (Reines & Volonteri 2015). However, the properties of dwarf AGN hosts are different from those of massive galaxy AGN hosts, which can complicate detecting AGN signatures. For example, most dwarfs are gas-rich (Kannappan 2004), strongly star-forming (Geha et al. 2012), and metal-poor (Tremonti et al. 2004). Supernova feedback preferentially expels metal-rich gas (Mac Low & Ferrara 1999), while accretion of low metallicity, intergalactic medium gas drives star formation (Dekel & Birnboim 2006).
Unlike in the unified AGN model, IMBHs in dwarfs often wander within 1 kpc of the center due to the dwarfs' weak gravitational potential (Reines et al. 2020, Bellovary et al. 2021. Therefore, it is unclear whether the AGN and stellar radiation fields strike the same gas clouds, or spatially separated gas clouds, as the IMBH relocates or settles down in a particular location. This uncertainty calls into question applying the centralized geometric model for massive AGN to dwarf AGN. Another issue for dwarf AGN, unlike massive AGN, is that X-ray observations that could provide valuable constraints on the AGN spectral energy distribution (SED) remain rare for M BH ≈ 10 5 M , and absent for M BH ≈ 10 3 M (Desroches et al. 2009). The few sources luminous enough to observe are likely outliers having fortuitous conditions to enable their detection (e.g. Godet et al. 2012). This situation introduces doubts about whether physical quantities controlling the shape of the SED scale down to the lowest black hole masses (Arcodia et al. 2020).
Photoionization models can provide the missing ingredient for IMBH detection by systematically accounting for uncertainties pertaining to the gaseous geometry and properties of low-mass AGN. The left panels of Figure 1 illustrate open, or plane-parallel geometries, where the covering factor is small, as typically assumed when modeling AGN (Elvis 2000;Feltre et al. 2016) and starforming regions like the Orion blister (Ferland 2001). The right panels of Figure 1 illustrate closed, or spherical geometries, where the covering factor is close to unity, as typically assumed when modeling obscured AGN like ULIRGs (Abel et al. 2009) and star-forming regions like 30 Doradus (Pellegrini et al. 2011).
Similarly, photoionization models have used two methods to take into account simultaneous AGN and stellar excitation, but rarely with justification. Figure 1, top panels, illustrates one approach where an AGN and starlight strike the same cloud and thus the SEDs from each source are mixed a priori (e.g. Abel et al. 2009), which we label as coincident mixing. The bottom panels of Figure 1 show another approach where the excitation sources illuminate spatially separated clouds and thus the models are mixed a posteriori (e.g. Kewley et al. 2013, Meléndez et al. 2014, Richardson et al. 2014), which we label as noncoincident mixing. Since dwarfs hosting IMBHs display a variety of morphologies (Kimbrell et al. 2021), it is unclear whether a particular gas geometry or mixing methodology would generically apply. Therefore, assessing all possibilities in Figure 1 is paramount.
Optical spectroscopy enables the categorization of emission line galaxies as AGN, star-forming, or a mixture of the two using [O III] λ5007/Hβ against [N II] λ6584/Hα (i.e., the BPT diagram), [S II] λ6720/Hα, and [O I] λ6300/Hα to form diagnostic diagrams (Baldwin et al. 1981, Veilleux & Osterbrock 1987. Previous photoionization modeling including AGN suggests that optical lines might grow too faint for detection for black holes outside 10 6 − 10 9 M , thus skewing M BH distributions (Cann et al. 2019, Bhat et al. 2020; however, none of this modeling accounts for the multiple geometrical configurations indicated in Figure 1, for the presence of stellar excitation, or for the uncertainty in the shape of the IMBH SED. Indeed, optical observations may contradict the theoretical impossibility of detection of IMBH AGN. Broadline selected dwarf AGN with M BH ≈ 10 5 M can be sometimes optically classified as AGN using the BPT diagram (Barth et al. 2004), indicating high [O III]/Hβ. Reines et al. (2020) used radio interferometry to identify dwarf AGN with optical star-forming galaxies clas- and two different mixing methodologies (coincident and non-coincident). Dwarf AGN could presumably fall into any of these four categories and therefore photoionization models must account for this systematic uncertainty. Note that the rays of light (arrows) only show one possible path to the observer and do not represent all possible components of the observed spectrum.
sifications, and used the M * − M BH relation to deduce M BH ∼ 10 4.1 − 10 5.8 M , albeit this relation shows up to 1.0 dex scatter for M * < 10 9 M . Optically classified star-forming galaxies might contain a treasure trove of additional hidden IMBHs that evade detection on account of the BPT diagram preferentially identifying high metallicity AGN (Polimera et al., submitted; hereafter P21). However, active black holes in the 10 3 M regime still remain undetected.
IR spectroscopy offers a better opportunity to reliably detect IMBHs as it possesses a wealth of high ionization lines insensitive to gas metallicity and dust extinction. The Spitzer era revealed the potential of the mid-IR to separate starbursts from AGN, leading to the development of several diagnostic diagrams involving [Ne V] and [O IV] emission lines and polycyclic aromatic hydrocarbon (PAH) features (Dale et al. 2006). The James Webb Space Telescope (JWST), spanning 0.6−28.3 µm, has the potential to revolutionize the search for black holes at the low-mass end of the IMBH distribution. Indeed, recent work has shown the spectral range of JWST can uncover AGN eluding detection from optical spectroscopy (Satyapal et al. 2020), although the uncertainties stemming the from IMBH SED and the configurations presented in Figure 1 remain unexplored.
In this paper, we fill in the gap in photoionization modeling to account for the uncertainty in gaseous geometry, mixed excitation, and AGN SED shape. We make emission line predictions that will be valuable for searching for 10 3 − 10 5 M black holes with optical spectroscopy and future JWST observations, while freely providing our simulation suite to the community 1 .

Incident Radiation Field
Figure 2 displays the three different models for the AGN SED that we have explored, assuming M BH = 10 3 − 10 5 M : "disk-plaw," "Cloudy," and "qsosed." The "disk-plaw" SED combines the diskbb accretion disk model (Mitsuda et al. 1984) with a power law (Γ = 2.1), normalized to give α ox = 1.41 (Grupe et al. 2010) where The inner diskbb temperature is calculated using Peterson (1997), where we assumeṁ/ṁ Edd = 0.1 and R = 3R S . The result is a piecewise, physically motivated SED for photoionization modeling (Cann et al. 2019, Bhat et al. 2020, but lacks physical self-consistency, such as the accretion disk radiation generating seed photons for the hard X-ray component. The "Cloudy" SED is the Cloudy (Ferland et al. 2017) default AGN SED, which is an empirical model that assumes the observed SED is the same as the continuum seen by nebular clouds. The functional form of the SED is given by, where T peak = 0.77 T in (Mitsuda et al. 1984), α UV is the low energy slope in the UV, α X is X-ray slope, and a is a constant adjusted to satisfy α ox . All of the spectral indices are taken from the median of the extinction corrected BLS1s sample in Grupe et al. (2010). The "qsosed" uses the agnsed model (Kubota & Done 2018), while fixing most parameters to their "typ-10 16 3 × 10 15 4 × 10 15 6 × 10 15 2 × 10 16 [Hz] 10 1 10 2 2 × 10 1 3 × 10 1 4 × 10 1 6 × 10 1 h [eV]  ical" values and assumingṁ/ṁ Edd = 0.1, a * = 0, and i = 45 • . This is a physically self-consistent SED for photoionization modeling (Panda et al. 2019, Vasiliev et al. 2020, Sarkar et al. 2021) appropriate for sub-Eddington accretion in IMBHs. Figure 2 shows profound differences between the three SED models. In part, the differences are due to the selfconsistent physics in the qsosed that is not featured in the other two models. The qsosed selects the inner accretion disk radius needed to power the X-ray emission and then truncates the disk at that point, rather than arbitrarily selecting R in = 6GM/c 2 . This results in the qsosed SED for M BH = 10 3 M having approximately the same peak energy as the other two SEDs for M BH = 10 5 M , and also a more appreciable hard X-ray component. Throughout the remainder of the paper, we only use the two SEDs that represent extremes: disk-plaw and qsosed.
To model the starburst continuum, we use the binary stellar population synthesis (SPS) code BPASS v2.0 (Stanway et al. 2016) to predict the spectrum emitted from stars subject to binary evolution. Including binaries is essential for a realistic treatment of the Wolf-Rayet (WR) phase, which can result from mergers and envelope removal (D'Agostino et al. 2019 to the age where the ionizing continuum flux for He + and O ++ (>54 eV) reaches a maximum. After continuous star formation for 250 Myr, the ionizing continuum ceases to evolve and this provides a "general" use stellar SED for photoionization modeling. We include 11 different metallicities calibrated to Z = 0.02 spanning 0.05 Z -2.0 Z for all SFHs. Figure 3 displays stellar SEDs for each of the three SFHs at 0.4 Z .

Gaseous Cloud
Following Richardson et al. (2019), we select a hydrogen density of log n H = 2.0 [cm −3 ] at the illuminated face. After selecting n H , the ionization parameter is then given by, where φ H is the hydrogen ionizing flux. While most emission line galaxies indicate log U = -3.5 --2.0, an even lower limit is needed to explain the lowest ionization dwarf galaxies and log U > -1.5 is needed to explain local blue compact dwarf galaxies (Stasińska et al. 2015). Accordingly, we run simulations with ionization parameters from log U = -4.0 to log U = -0.5 in increments of ∆(log U ) = 0.25. We follow Abel et al. (2008) Table 1. The scaling of abundances with metallicity includes a detailed prescription for accounting for specific elemental variations due to differences in nucleosynthesis. Unfortunately, the abundance of oxygen relative to solar has become synonymous with metallicity. While metallicity is strictly defined as the mass fraction of metals, and oxygen makes the greatest contribution to metallicity, the two are not equivalent. To avoid this ambiguity, we refer to the scaling parameter ζ O as the metallicity (see Nicholls et al. 2017) where the solar metallicity is 12+log(O/H) = 8.76, corresponding to ζ O = 1.
We assume Orion grains and PAHs throughout the cloud as implemented in Baldwin et al. (1991) and Abel et al. (2008), respectively. The dust abundance is typically assumed to not depend on metallicity, while the gas (hydrogen) to dust ratio (G/D) varies as G/D ∝ Z −1 (Dwek 1998). Here, we adopt a more sophisticated broken power law, (5) which shows that low metallicity dwarfs deviate from a single power law relation (Rémy-Ruyer et al. 2014).
The Orion grain abundances and PAH abundances are scaled from their default abundances by the same factor to satisfy this relationship at a given Z.
Most photoionization modeling assumes a "standard" set of gas phase depletion factors derived from a variety of sources. In reality, however, depletion factors δ X depend on the undepleted reference abundance set being used and the strength of the depletion F * . We follow a self-consistent approach where we use our unique reference abundances with the methodology outlined in Jenkins (2009). The strength of the depletion is selected so that the depletion factor for iron is δ F e = −1.5 dex as in Thomas et al. (2018). This results in δ O = −0.11 dex, which is less than the commonly assumed Cloudy default of δ O = −0.22 dex, but assists in matching optical emission line diagnostics ) and matches the depletion derived from analyzing dust grain composition (Peimbert & Peimbert 2010), X-ray spectroscopy (Pinto et al. 2013), and α-element enhancement (Amayo et al. 2021). In reality, depletion factors should change with G/D at a fixed metallicity and as a function of metallicity (Peimbert & Peimbert 2010, De Cia et al. 2016; however, this is rarely taken into account in photoionization modeling and beyond the scope of this work. The complete set of adopted depletion factors are listed in Table 1 and we elaborate on our methodology in Appendix A. Our final gas phase abundances differ from the BPASS stellar abundances, the effects of which have been investigated in Grasha et al. (2021), but are beyond the scope of this work. All in all, our model suite consists of > 6.33 × 10 4 simulations.

LINE RATIO SENSITIVITY ANALYSIS
We seek to assess the sensitivity of emission line ratios to f AGN , AGN SED shape, mixing methodology, and geometry. In particular, we select emission line diagnostics that are detectable in purely star forming galaxies, enabling investigation of the connection between star formation and AGN in dwarfs, which is not possible for all emission lines. For example, the presence of [Ne V] 14.3µm or [Ne VI] 7.65 µm alone signals AGN activity, but our simulations show that [Ne V] is unlikely to be detectable in most local dwarfs where U and f AGN are small, thus making it an unreliable tracer of AGN fraction for all galaxy masses.
Similarly, coronal lines from highly ionized states (e.g., Si VI, Fe XIII) have been used to identify AGN (Cann et al. 2018, Bohn et al. 2021, Kimbro et al. 2021). However, several limitations exist to this approach: (1) only about half of AGN actually show coronal line(s), regardless of instrumental line sensitivity (Riffel et al. 2006); (2) highly ionized states depend on physical conditions with high U ∼ −2.0 and high f AGN > 0.64, which are not characteristic of typical dwarfs; (3) if present, the coronal line region lies between the BLR and NLR, which implies dust sublimation is important, a process that presents a problem for self-consistent photoionization modeling (Mazzalay et al. 2010, Adhikari et al. 2016); (4) coronal line emission typically originates from metals that become heavily depleted in forming dust grains (e.g., Si, Ca, Fe), and therefore small changes to F * yield large differences in these abundances (De Cia et al. 2016). All together, these limitations suggest that using coronal lines will be subject to strong selection bias, and therefore alone cannot provide a complete picture of IMBH activity in dwarf AGN.
With these caveats in mind, we have determined three different mass ranges to evaluate the AGN emission line diagnostics within the wavelength ranges of SDSS and JWST: lower mass dwarfs, intermediate mass dwarfs, and massive galaxies.
• For lower mass dwarfs, we use the U −Z correlation presented in Kashino & Inoue (2019), which leads to values of log U ≈ −2.0, Z/Z ≈ 0.15.
• For intermediate mass dwarfs, we analyzed SDSS strong emission line measurements (Tremonti et al. 2004) as processed by P21 for the z∼0 dwarf dominated RESOLVE Survey (Kannappan & Wei 2008) with the Bayesian analysis code NebulaBayes (Thomas et al. 2018), but using our simulation suite. We determine median values of log U ≈ −3.25, Z/Z ≈ 0.4 for intermediate mass dwarfs (P21), which can be thought of as a 'fiducial' set of parameters for dwarf galaxies.
• For massive galaxies, we adopt log U ≈ −2.0, Z/Z ≈ 1.0 following the parameters used in (Abel et al. 2009) to model AGN in massive galaxies such as ULIRGs.
For the rest of §3, we use the stellar SED resulting from continuous SFH at 250 Myr as it represents a "general purpose" ionizing continuum not purposefully constructed for maximizing any particular emission line ratio.  The photoionization cross section of hydrogen quickly decreases as photon energy increases (Osterbrock & Ferland 2006). As a result, X-rays due to an AGN penetrate into neutral gas, which causes collisions that excite neutral species. Since [O I] λ6300 partially originates from neutral gas, [O I]/Hα traces f AGN rather well in dwarf galaxy mass regimes in contrast to the other optical diagnostics. This makes it a reliable choice for detecting AGN activity in dwarfs and constraining f AGN . Nonionizing photons can also pump electrons into high energy excited states, which can cause [O I] λ6300 emission when the electrons cascade to lower energy levels. This continuum fluorescence process is rather inefficient for O I (Bautista 1999), leaving collision excitation as the dominant process, and thus preserving the utility of [O I] λ6300 as an AGN diagnostic.

Optical
The AGN SED and mixing methodology can create substantial variation in all diagnostics, although this is mitigated for emission line ratios created under certain physical conditions. For example, the middle rows of from the metallicity sensitive BPT diagram, corroborates that dwarf AGN have inconsistent optical classifications (Reines et al. 2020, P21) for black hole masses in the range 10 3 − 10 5 M .
Our results also stand in contrast to other work (Cann et al. 2019), suggesting that AGN host galaxies with

Mid-IR
The most valuable AGN diagnostics for JWST come from the MIRI instrument (4.9 -28.3 µm), as opposed to NIRspec (0.6 -5.3 µm), based on the availability of either emission lines originating from the same element with different ionization potentials, or high ionization lines along with a recombination line for comparison. Although we are focused on detecting active IMBHs in relatively local galaxies, it is worth noting that galaxies at z > 0.3 will have optical diagnostics that fall within NIRspec, thereby providing additional constraints.  become less pronounced as AGN fraction increases; in contrast to the optical, IR line ratios generally show less sensitivity to geometry. Weaver et al. (2010)   In general, line ratios featuring [O IV] and [Ar II] show less sensitivity to geometry than optical diagnostics, display more variation with AGN fraction than other mid-IR line ratios, and retain their diagnostic potential over a wider range of galaxy masses and BH masses. The mid-IR diagnostics provide greater constraints on f AGN for lower levels of AGN activity expected for dwarf AGN with IMBHs, compared to optical diagnostics which are more sensitive to higher levels of AGN activity. These results highlight the benefit of using multiple line ratios to constrain f AGN with diagnostic diagrams.

DIAGNOSTICS DIAGRAMS
In light of the wide range of uncertainties previously mentioned, we seek to find diagnostic diagrams suit- able for detecting AGN excitation, and ideally measuring f AGN , for a given galaxy. To compare our models to galaxies with spectroscopic observations, we select a subset of models (−3.5 ≤ log U ≤ −1.5) with an instantaneous SFH, since the harder stellar continua provide a limiting case for AGN activity.

Optical
We use a 20 Myr SSP for the optical diagnostics since this SFH maximizes the [O III]/Hβ line ratio used in many diagnostic diagrams that separate AGN and star forming galaxies. For the observational sample, we choose dwarf galaxies in the MPA-JHU sample from P21. Following P21, we classify dwarf AGN as any dwarf galaxies not classified as definitely star forming (e.g., traditional AGN, LINERs, etc.) The three diagrams from Veilleux & Osterbrock (1987) are by far the most commonly used optical diagnostics for separating AGN and star forming galaxies, so we focus our analysis on them.  Figure 1. This large spread in values relates to the current difficulty of detecting IMBHs at this mass, but also provides promise that optical detection is indeed possible. Figure 8 shows that models in the range 0 ≤ f AGN ≤ 0.16 occupy similar regions of each diagram. In particular, the BPT diagram (top row) registers observed dwarf AGN as consistent with pure star formation. Conversely, theoretical models with pure star formation can cross over into the AGN region of the diagram. Despite the differences in photoionization modeling, this is also seen in Xiao et al. (2018) on account of using the BPASS SEDs that feature substantially harder continua than many other SPS codes.
The [S II]/Hα diagram indicates that pure star formation can reproduce the line ratios of all the dwarf AGN shown, but unlike for the BPT diagram, for the [S II]/Hα diagram this result only occurs with the par- Figure 8. Optical diagnostic diagrams that separate AGN activity from starburst activity assuming MBH = 10 3 M (left column) and MBH = 10 5 M (right column). The selected age of the binary stellar population maximizes the contribution of WR stars and thereby sets a lower limit for AGN activity. Pure starburst models (0% AGN) are displayed as blue stars. The models for a given fAGN approximately span the area of each hatched shape. The demarcations in each diagram are taken from Kewley et al. (2001), Kauffmann et al. (2003), and Kewley et al. (2006). The majority of observed dwarf AGN lie in the star forming wing of the BPT diagram (top panels), which models with and without AGN can reproduce.
ticularly hard stellar SED that we have selected. Thus, knowing the SFHs for a given sample of dwarfs could provide justification for using [S II]/Hα to identify AGN. In contrast, the BPT diagram is more generally a poor diagnostic for identifying dwarf AGN, regardless of SFH, as discussed in more detail in P21.
The last optical diagnostic diagram shows that some dwarf AGN require models with an AGN component to reproduce the observed [O I]/Hα, making this diagram more reliable in separating star formation and AGN in dwarfs (Reines et al. 2020, P21). It is worth noting that more dwarf AGN likely require an AGN component than shown here, on account of the hard stellar SED used in this analysis. Additionally, for [O I]/Hα there is some modest separation between models with different f AGN , especially with M BH = 10 3 M , enabling a better assessment of AGN activity as also shown in Figure 4.

Mid-IR
For our observational sample, we choose dwarf SF (Cormier et al. 2015) and dwarf AGN (Hood et al. 2017) observations requiring that the AGN satisfy [Ne V]/[Ne II]>0.1 (Inami et al. 2013). It is important to note that Pfα, featured in Figure 6, is not available in statistical samples due to the low sensitivity and resolution of Spitzer and ISO. Similarly, since [Ar II] 6.98µm is unavailable for samples including dwarfs, we use massive AGN (Sturm et al. 2002) for diagnostic diagrams with this line.
We have tested all possible emission line ratio combinations from the emission lines listed in §3.1.2. Figure 9 shows the results of this analysis in three diagnostic diagrams featuring emission line ratios also presented in Weaver et al. (2010), Inami et al. (2013), andHao et al. (2009).
Our new demarcations separate three distinct regions: the SF region contains only pure starburst galaxies; the SF and/or AGN region contains a mixture of pure starbursts and AGN; the AGN region contains only AGN. The divisions between SF and SF and/or AGN regions are given by,  (Hao et al. 2009) and [S IV]/[Ar II] in the abscissa, which creates even more separation between starbursts and AGN, and also between f AGN values. As mentioned above, the [Ar II] observations are poor since Spitzer was only capable of detecting wavelengths <10 microns in low resolution mode (R ∼ 130). In contrast, MIRI features the highest resolution (R ∼ 3100) and line sensitivity available at these wavelengths. Thus, the right panels in Figure 9 present a promising new diagnostic for JWST observations to constrain AGN activity.

DISCUSSION AND CONCLUSIONS
In this paper, we have explored several key uncertainties that will be critical to assess in the upcoming era of detecting IMBHs: the shape of the AGN SED, the surrounding gas cloud geometry, and the manner of mixing radiation from AGN and stars. The diversity of models for the AGN SED in Figure 2 emphasizes the importance of next generation X-ray facilities such as Lynx and Athena that will enable proper constraints on low luminosity sources. We show that all of the listed uncertainties have profound effects on emission line diagnostics in the optical and mid-IR. This has implications for statistical samples of galaxies that are analyzed with photoionization models that implicitly assume a particular geometrical configuration or mixing methodology in Figure 1. Such assumptions can lead to powerful selection effects (Ferguson et al. 1997, Meskhidze & Richardson 2017. These assumptions could limit the applicability of simple relationships that scale physical properties (e.g., f AGN , M BH ) with emission lines ratios for large samples of galaxies.
Even though degeneracies abound in the optical, we find AGN diagnostics remain detectable in most situations across a range of galaxy masses. This result offers promise that surveys like SDSS have already picked up signatures of AGN with M BH = 10 3 M . Unfortunately, the gold standard for optical AGN classification, the BPT diagram, poorly traces AGN activity in dwarf galaxies. As we argue in P21, this is due to [N II]/Hα Figure 9. Mid-IR diagnostic diagrams that separate AGN activity from pure starburst activity in the same format as Figure 8. Pure starburst models (0% AGN) are displayed as blue stars. The models for a given fAGN approximately span the area of each hatched shape. The SF, SF and/or AGN, and AGN regions represent pure starburst galaxies, a mix of pure starburst galaxies and AGN, and only AGN, respectively (eqns 6, 7, 8). JWST will enable emission lines normally missing from Spitzer spectra to serve as additional constraints, like [Ar II] 6.99 µm (left panels).
being metallicity sensitive and f AGN insensitive. It also reflects the sensitivity of [O III]/Hβ to the uncertainties explored in this paper. Actually, as seen in Reines et al. 2020 and P21, the often forgotten [O I]/Hα diagnostic is a stronger metric for an active AGN in dwarfs given its sensitivity to f AGN and relative insensitivity to physical uncertainties. It therefore serves as the best optical emission line ratio for finding dwarf AGN, but the relative weakness of the [O I] line does limit the sample that can be tested (P21).
Relatively high-mass dwarfs (M * 10 9.3−9.5 M ) classified as LINERs present an opportunity to find the coveted 10 3 M black holes, which would fill in a major gap in IMBH empirical relationships and help constrain the occupation fraction in dwarfs. If AGN in this BH mass regime are confirmed by other methods (e.g., variability, line broadening, X-rays), but still fail to show any optical signatures, this non-detection would suggest that the column density of obscuring gas and the covering factor around the source could be factors in the elusive nature of these AGN. The column densities and covering factors of dwarf AGN remain relatively unexplored. Purely star forming dwarf galaxies show a decrease in covering factor as metallicity decreases, albeit with a large scatter in the relation (Cormier et al. 2019). Conversely, massive galaxies show a decrease in obscured AGN (large covering factor) with increasing luminosity (Sazonov et al. 2015, Georgakakis et al. 2017, and therefore increasing metallicity (Lamareille et al. 2004). It is possible that dwarf AGN could provide the missing link between these two opposing trends if low mass IMBHs continue to remain elusive in the optical.
In the mid-IR, we have shown that [Ar II] 6.98µm / Pfα and [O IV] 25.9µm/Pfα together have the potential to constrain AGN activity over a range of galaxy masses while minimizing sensitivity to physical uncertainties ( Figure 6) (Figure 9). These new diagrams are capable of separating starburst activity from AGN activity down to a 4% AGN fraction in dwarfs, which are ubiquitously star forming. Unlike the BPT diagram, these proposed diagrams maintain their diagnostic value over a wide range of ionization and metallicity.
We emphasize that the provided demarcations serve as theoretical boundaries for classifying the excitation mechanism in dwarf galaxies. Observations of optically classified starbursts can occasionally cross over into the AGN region of these diagrams due to strong mid-IR [O IV] emission (e.g., IZw18, Lebouteiller et al. 2017), although soft X-ray and hard X-ray emission is typically detected in these galaxies.
Current SPS models are unable to produce enough photons >54 eV to account for this crossover. High mass X-ray binaries (HMXBs) can generate the hard photons needed to create AGN-like emission line ratios assuming a model with high U , high L X /SFR, and a sufficiently hard SED (Simmonds et al. 2021). It is unlikely that HMXBs can uniformly account for dwarfs exhibiting emission line ratios characteristic of AGN given that the shape of the HMXB SED remains highly uncertain and most dwarf AGN do not present these extreme conditions (P21). Additionally, a more generic treatment of the HMXB SED has shown that HMXBs are an inefficient means of producing photons >54 eV (Senchyna et al. 2020). Therefore, unless more realistic treatment of the WR phase alters the result of our predictions in the future, our current simulations suggest these galaxies have active AGN. The SF and/or AGN region in each diagram can contain optically classified starbursts (see Figure 9), and this result is reproduced by a photoionization model with multiple gas clouds (Meléndez et al. 2014, Richardson et al. 2016. IFU observations should help clarify whether such a model needs to be invoked. To fully address modeling degeneracies in both the optical and mid-IR, Bayesian analysis will be a valuable tool for extracting meaningful properties over a range of conditions in future studies. Our simulation suite incorporates a finer spacing of metallicities, ideal for accurate Bayesian analysis at < 0.4 Z (Richardson et al. 2019), than other SPS models (e.g., Starburst99), while robustly accounting for D/G and elemental depletion. In addition, the mixing methodologies and geometries unique to our models could be used in a multicomponent ISM analysis, similar to work on the multiphase ISM in dwarfs (Cormier et al. 2019). Together with the AGN diagnostics and excitation diagrams presented in this paper, such future applications of our models will prove useful for addressing more complex topologies associated with mixed AGN/SF excitation in dwarfs in the era of JWST.

ACKNOWLEDGMENTS
CR gratefully acknowledges the support of the Elon University FR&D committee and the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. This work used the XSEDE resource Comet at the San Diego Supercomputing Center through allocation TG-AST140040. MP gratefully acknowledges the support of the 2020 Hamilton Award from the UNC Department of Physics and Astronomy. SK and MP acknowledge support from NSF AST-2007351. JB acknowledges grant support from NSF AST-1812642 and CUNY JFRASE. We thank Nick Abel, Jenna Cann, Chris Done, Gary Ferland, Ed Jenkins, and Shobita Satyapal for helpful discussions that improved the quality of this paper.
Despite the dependence on a reference set of abundances, most photoionization modeling does not properly take this factor into account as shown below. Different lines of sight yield different depletion factors for a given element, allowing one to define a depletion strength, F * , that accounts for this variation. A linear fit to the logarithm of the depletion factor gives the following form, where A X is the slope, B X is the vertical offset, and z X accounts for the errors in the observations. We assume that the non-refractory elements He, Ne, S, and Ar do not become depleted in the ISM. However, sulphur depletion remains a subject of debate ( We use the fits for the Galaxy provided in J09 for C, N, O, Mg, Si, P, Cl, Ti, Cr, Mn, Fe, Ni, Cu, and Zn. We make this choice for two reasons. First, J09 provides the most complete sample, while extragalactic fits are mostly limited to heavier elements that have less impact on emission line predictions. Second, the number of observations included in the J09 analysis makes the fit more reliable. While we assume that depletion occurs as observed in the Milky Way, deviations in the depletion patterns for heavy elements are known to be present for other galaxies. For example, Mg, Ti, and Mn show noticeable deviations in the Small Magellanic Cloud (Jenkins & Wallerstein 2017), and while the depletion pattern is mostly Galactic in the Large Magellanic Cloud (Tchernyshyov et al. 2015), Si does show a noticeable deviation (Roman-Duval et al. 2019).
For the rest of the elements, we have compiled a list of depletion factors towards the highly depleted star ζ Oph to serve as the value at F * = 1.0 in accordance with the method used by J09. Fluorine is the only exception where we use the most depleted source in Snow et al. (2007), since the value towards ζ Oph is unavailable. To determine the parameters A X and B X for these elements, we rely on the general trend that elements with larger condensation temperatures show greater depletions and steeper slopes. The following elements have similar depletion trends to their J09 counterparts: Li ⇔ Cr; Be ⇔ Mg; B ⇔ Zn, F ⇔ Cl; Na, K ⇔ Cu, Zn; Al, Ca, V ⇔ Ti; Sc ⇔ P, Cu; Co ⇔ Ni. From these analogies, we can assume that all of the elements in J09 have δ X = 0 at F * ≈ −0.5 except for δ Fluorine , which reaches zero at F * = 0.2. Table 2 lists the depletion factors for all elements for F * = 1.0, the linear fit parameters need to recreate the depletion pattern, and the reference for each depletion factor given. Note that depletion factors have been rescaled according to our reference abundances instead of using the abundance assumed in each cited source. Displaying the depletion factors for F * = 1.0 makes it clear which elements are weakly depleted in the ISM as opposed to not depleted at all. However, for our photoionization modeling, we adjusted F * so that δ Fe = −1.5 (see Table 1), which we justified in §2.2. The methodology we have used here is similar to what resulted in the depletion factor sets included with the photoionization code Mappings V (Sutherland & Dopita 2017) with two major differences. First, unlike the values included with Mappings, our depletion factors are self-consistently scaled with our reference abundances. Second, we do not deplete nitrogen for any value of F * given the lack of evidence that nitrogen gets locked up in grains. This is particularly important for any analysis that strongly relies upon nitrogen emission lines for inferring metallicity values.
The depletion factors we use represent a small step toward accounting for the different depletion patterns in gas phase abundances, much in the same way Nicholls et al. (2017) accounts for nucleosynthesis patterns for elements as a function of metallicity. However, obvious limitations exist since galaxies can show more negative F * than we account for here (Jenkins & Wallerstein 2017). In this section, we present the line ratio sensitivity for additional AGN diagnostics. Figures 10 and 11 display diagnostics in the optical for M BH = 10 3 M and M BH = 10 5 M , respectively, and similarly in the mid-IR for Figures 12 and 13. The columns of each figure represent a given galaxy mass, while the rows represent a given line ratio. As with the diagnostics presented in the main body of the paper, the uncertainty introduced from the AGN SED shape, cloud geometry, and mixing methodology, decreases at higher black hole mass.
In the optical, most of the emission line ratios poorly trace f AGN . It is noteworthy that [N II]/Hα remains relatively insensitive to AGN activity in most conditions, which combined with the line ratio's metallicity sensitivity makes it an overall poor diagnostic for dwarf AGN (Figures 6 and 7) across the black hole masses and galaxy masses we have considered.