# Nebular emission from composite star-forming galaxies - I. A novel modelling approach

Christophe Morisset, Stéphane Charlot, Sebastián F. Sánchez, Carlos

Espinosa-Ponce, Eric Barat, Thomas Dautremer

► To cite this version:

Christophe Morisset, Stéphane Charlot, Sebastián F. Sánchez, Carlos Espinosa-Ponce, Eric Barat, et al.. Nebular emission from composite star-forming galaxies - I. A novel modelling approach. Monthly Notices of the Royal Astronomical Society, 2025, 538 (3), pp.1884-1905. 10.1093/mnras/staf143. insu-05101667

HAL Id: insu-05101667

https://insu.hal.science/insu-05101667v1

Submitted on 8 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Nebular emission from composite star-forming galaxies – I. A novel modelling approach

Christophe Morisset<sup>D 1,2★</sup> Stéphane Charlot<sup>D 3</sup>, Sebastián F. Sánchez<sup>D 1</sup>, Carlos Espinosa-Ponce<sup>D 4</sup>, Eric Barat<sup>D 5</sup> and Thomas Dautremer<sup>D 5</sup>

Accepted 2025 January 8. Received 2025 January 8; in original form 2024 October 8

## ABSTRACT

We introduce a novel approach to modelling the nebular emission from star-forming galaxies by combining the contributions from many H II regions incorporating loose trends in physical properties, random dust attenuation, a predefined Hα luminosity function, and a diffuse ionized-gas component. Using a machine-learning-based regression artificial neural network trained on a grid of models generated by the photoionization code CLOUDY, we efficiently predict emission-line properties of individual H II regions over a wide range of physical conditions. We generate 250 000 synthetic star-forming galaxies composed of up to 3000 H II regions and explore how variations in parameters affect their integrated emission-line properties. Our results highlight systematic biases in oxygen-abundance estimates derived using traditional methods, emphasizing the importance of accounting for the composite nature of star-forming galaxies when interpreting integrated nebular emission. Future work will leverage this approach to explore in detail its impact on parameter estimates of star-forming galaxies.

**Key words:** software: machine learning – software: simulations – HII regions – galaxies: abundances.

## 1 INTRODUCTION

Emission lines in the integrated spectra of star-forming galaxies provide valuable clues about the physical conditions in the interstellar medium (ISM) heated by young massive stars, such as gas metallicity, electron density and temperature, ionization conditions, and dust content. In practice, interstellar parameters, particularly gas-phase metallicities, of spatially unresolved galaxies have been determined using various emission lines at ultraviolet, optical, and infrared wavelengths using different techniques. Common approaches include the ‘direct method’ based on electron-temperature estimates from auroral lines (e.g. Aller & Liller 1959; Peimbert 1967; Peimbert & Costero 1969; Pérez-Montero 2017), ‘empirical’ calibrations linking auroral-line metallicities to strong-line flux ratios (e.g. Jensen, Strom & Strom 1976; Alloin et al. 1979; Pagel et al. 1979; Kewley & Ellison 2008; Marino et al. 2013), and fully theoretical approaches using photoionization models (e.g. Charlot & Longhetti 2001; Kewley & Dopita 2002; Morisset et al. 2016; Byler et al. 2018; Pérez-Montero et al. 2021, see Fernández et al. 2023 for a comparison between the direct method and a model-based method). Each method has its own limitations and can lead to significant discrepancies between the derived measurements (see the review by Kewley, Nicholls & Sutherland 2019).

Another complication in estimating interstellar parameters from the integrated nebular emission of star-forming galaxies is that these galaxies typically harbour a population of H II regions with different temperatures, ionization parameters, and metallicities. Yet, all the above standard methods are generally applied under the assumption that the ISM conditions in star-forming galaxies can be described by global parameters at the galaxy scale. Kobulnicky, Kennicutt & Pizagno (1999) showed that, while analyses of global emission-line spectra using standard methods provide fairly reliable constraints on chemical properties, they introduce small systematic errors due to varying gas conditions. The pioneering analysis by Kobulnicky et al. (1999) was based on a limited set of six combinations of two spectra from H II regions of fixed metallicity sampling three values of electronic temperature and two values of $O^+/O^{2+}$. Mast et al. (2014) explore how spatial resolution affects observed line ratios and derived oxygen abundances, and how integrating the emission of entire galaxies biases spectral analyses. Also relevant is the introduction by Ferguson et al. (1997), to model the spectra of Seyfert galaxies, of the concept of ‘locally optimally emitted cloud’, where a distribution of clouds with a wide range of densities and distances to the ionizing source is used to reproduce the narrow-line emission. This concept has been used by Richardson et al. (2016) and Lebouteiller et al. (2024) in the case of star-forming galaxies. Recently, Cameron, Katz & Rey (2023) used a radiation-hydrodynamics simulation of an isolated dwarf galaxy to show that temperature fluctuations in the ionized gas throughout the galaxy lead to an underestimation of

© 2025 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1">¹Universidad Nacional Autónoma de México, Instituto de Astronomía, AP 106, 22800 Ensenada, BC, México</footnote>

<footnote n="2">²Universidad Nacional Autónoma de México, Instituto de Ciencias Físicas, Av. Universidad s/n, 62210 Cuernavaca, Mor., México</footnote>

<footnote n="3">³Sorbonne Université, CNRS, UMR7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="4">⁴Universidad Nacional Autónoma de México, Instituto de Astronomía, A. P. 70-264, C.P. 04510 México, D.F., México</footnote>

<footnote n="5">⁵Université Paris-Saclay, CEA, LIST, F-91120 Palaiseau, France</footnote>

<footnote>*E-mail: morisset@astro.unam.mx, chris.morisset@gmail.com</footnote>
the true metallicity in estimates of O/H from global emission using the direct method. This type of simulation is, however, very time-consuming and cannot be easily performed for hundreds of thousands of galaxies. Yet, our understanding of the populations of H II regions in nearby star-forming galaxies has rapidly advanced in recent years, largely thanks to integral-field spectroscopic observations (e.g. Marino et al. 2013; Rousseau-Nepton et al. 2018; Santoro et al. 2022).

In this paper, we introduce a novel approach to model the nebular emission from star-forming galaxies hosting large populations of H II regions spanning wide ranges of physical properties. Our approach relies on the efficiency of machine learning (ML) to interpolate an unlimited number of models from a carefully sampled reference grid encompassing extensive ranges of H II-region parameters.

A wide range of ML-based algorithms can be used to perform interpolation in high-dimensional spaces, spanning from basic linear interpolation to more sophisticated techniques. For example, Galliano (2018) used linear interpolation to predict dust models, while Ramambason et al. (2022) applied a nearest-neighbour approach to interpolate within a grid of half a million models. Wu & Boada (2019) used radial-basis functions to interpolate within a grid of 1367 photodissociation-region models, and Smirnov-Pinchukov et al. (2022) employed a k-nearest neighbour regressor to interpolate within a grid of 1 million thermochemical protoplanetary-disc models. Furthermore, Bron et al. (2021) used random-forest regressors to predict observables in a chemistry-model grid, whereas Maltsev et al. (2024) leveraged Gaussian-process regression to emulate stellar tracks in a five-dimensional parameter space. Artificial neural networks (ANNs) have also been explored, with Ho (2019) employing them to interpolate within a grid of emission-line observations to derive oxygen abundances, Palud et al. (2023) investigating the chemical composition and emission of diffuse clouds and photodissociation regions, and Li et al. (2024) developing the Cue tool to predict stellar and nebular emission using ANNs trained on a set of CLOUDY<fnref n="1" /> models (Ferland et al. 2017). Autoencoder-based techniques (a subclass of ANNs) have also gained prominence, with Badenas-Agusti et al. (2024) applying them to interpolate within a grid of 25 000 stellar-atmosphere models, and Holdship et al. (2021) and Grassi et al. (2022) using autoencoders to emulate thermochemistry models and chemical networks, respectively.

In this work, we opt for a regression ANN with a simple architecture due to its exceptional speed during the prediction phase, which is crucial for our large-scale production of H II-region models to be incorporated in synthetic galaxy-emission models.

The method presented in this paper involves developing, training, and testing a regression ANN to predict model properties in record time for any set of input parameters. ANNs are a class of algorithms with architectures loosely inspired by the biological neural networks found in animal brains. Their adaptability and ability to model non-linear relations make them highly effective for a broad range of tasks, such as classification, regression, and clustering, which has contributed to their widespread adoption in

scientific fields, in particular in astrophysics. In this study, we use ANNs as regressors: predicting one or more floating-point outputs from a set of floating-point inputs. The ANN model is trained by optimizing its parameters (weights and biases of the individual neuron-like functions) through an iterative process, which minimizes the difference between the predicted outputs and the actual outputs in the training set. For more details on the use of ML in astrophysics, and a comprehensive description of ANN characteristics, see Baron (2019).

We use this method to compute the emission-line properties of 250 000 synthetic star-forming galaxies, obtained by combining the contributions from 100 to about 3000 H II regions per galaxy, incorporating loose trends in physical properties, random dust attenuation, a predefined Hα luminosity function, and a diffuse gas component ionized by Hα Low-Mass Evolved Stars (HOLMES; see Flores-Fajardo et al. 2011). Our method allows exploration of how different parameter distributions affect the integrated emission-line properties of composite galaxies. As an example of application, we study systematic biases in O/H estimates using the direct method, particularly at low and high metallicities. Our results highlight the importance of accounting for the composite nature of star-forming galaxies when interpreting integrated nebular emission, a topic we will explore further in the next paper in this series.

We present the reference grid of photoionization models of H II regions in Section 2 and the ANN and its performance in Section 3. In Section 4, we describe our approach to computing the nebular emission from synthetic star-forming galaxies. The properties of these galaxies are presented in Section 5, together with our study of the systematic biases implicit in O/H estimates using the direct method. Section 6 summarizes our conclusions.

## 2 AN ORIGINAL GRID OF NEBULAR-EMISSION MODELS OF H II REGIONS

The high dimensionality of the parameter space we wish to explore prevents us from adopting a regular Cartesian grid, the realization of which would require unmanageable resources. We therefore adopt a different approach, which consists in (i) running an irregular grid of photoionization models, where the different adjustable parameters can take random values drawn from selected distributions; (ii) using ML techniques to train, based on this grid, a model able to predict emission-line ratios of interest for given input values of the adjustable parameters; and (iii) generating in record time, with this ML model, grids of $10^8$ or more models for any (regular or not) distribution of adjustable parameters. The use of an irregularly sampled model grid for the training set has a number of advantages. It avoids oversampling of single parameter values with no influence on certain observables, minimizes potential discretization effects in interpolation, and allows the grid size to be increased more easily than with a Cartesian grid.

In this section, we describe our approach to producing a comprehensive grid of models. We start by defining the set of adjustable parameters in Section 2.1, with particular attention to chemical abundances and ionizing spectra. Then, in Section 2.2, we define the distributions of these parameters.

### 2.1 Irregular model grid in high-dimensional parameter space

We compute photoionization models of H II regions using the multipurpose photoionization code CLOUDY in its version 17.03.

<footnote n="1">¹CLOUDY is a versatile photoionization code designed to simulate the physical and chemical conditions in interstellar matter across of wide range of environments. It computes radiation transfer through a region photoionized by a user-defined spectral energy distribution, assuming thermal and ionization equilibrium. The code provides a variety of outputs, including emission-line intensities, ionic fractions, and electron temperature (see https://nublado.org for more details).</footnote>
They are defined, run, and stored under the reference 'HII_24' in 3MdB_17,² a database using the PYCLOUDY facility (Morisset 2014).

The adjustable parameters of the CLOUDY models considered in this work are

(i) the hydrogen density $n_H$;

(ii) the gas filling factor $\epsilon$;

(iii) the ionization parameter $U$: unless otherwise specified, $U$ refers here to the mean ionization parameter over the volume $V$ of the $H^+$ region, $U = \bar{U} = \frac{1}{V} \int_V U(r)4\pi r^2 dr$, with $U(r) = Q_0/(4\pi r^2 n_{Hc})$ at radius $r$. In this expression, $Q_0$ is the rate of ionizing photons, and $c$ is the speed of light. For a nebula with constant $n_{Hc}$ confined between $R_1$ and $R_2$ (as in the models considered here), $\bar{U} = 3U(R_2)(1 - R_1/R_2)/[1 - (R_1/R_2)^3]$;

(iv) the form factor $f_r$, defined as the ratio $f_r = R_{in}/R_{s}$ of the model inner radius to the Strömgren radius corresponding to $R_{in} = 0$. This parameter characterizes the shape of the nebula, transitioning from fully spherical for $f_r \ll 1$ to plane-parallel for $f_r \gg 1$ (see section 4.1 of Stasińska et al. 2015, where this parameter is called $f_s$). The relation between $\bar{U}$, $Q_0$, $n_H$, $\epsilon$, and $f_r$ for a pure-hydrogen nebula is $Q_0 = 4\pi c^3 \bar{U}^3 / \{3\alpha_B^2 n_H \epsilon^2 [(f_r^3 + 1)^{1/3} - f_r]^3\}$, where $\alpha_B$ is the case-B recombination coefficient;

(v) the chemical composition, i.e. the abundances relative to hydrogen of the 28 heavy elements handled by CLOUDY (see Section 2.1.1);

(vi) the parameter $\xi_d$, as defined by Gutkin, Charlot & Bruzual (2016), which represents the mass fraction of heavy elements depleted on to dust grains (see Section 2.1.1);

(vii) the parameter HβFrac: each individual CLOUDY model generates 10 entries in the data base, after being sliced at different optical depths corresponding to matter-bounded cases. The depth of the cut is parametrized by the quantity HβFrac, defined as the ratio of the Hβ luminosity of the truncated model to that of the radiation-bounded model. Thus, HβFrac < 1 corresponds to matter-bounded models and HβFrac = 1 to radiation-bounded models (see section 4.1 of Stasińska et al. 2015). The parameter HβFrac can be approximately related to $f_{esc}$, the fraction of ionizing photons escaping from the nebula, by the expression $f_{esc} = 1 - H\beta\text{Frac}$ (we also consider below, in Section 4.1, the 'picket-fence' model with gas-covering factor less than unity, where photons can escape through holes in the nebula);

(viii) the age $t_c$ and metallicity $Z$ of the ionizing star cluster, which determine the shape of the ionizing spectrum (see Section 2.1.2 for more details). In this work, we adopt the same metallicity for the stars as for the gas they ionize.

### 2.1.1 Abundances, depletion, and dust

We start by drawing the total (i.e. in both the gas and dust phases) oxygen abundance of a model, $(O/H)^{\text{tot}}$. The abundances of all other chemical elements are drawn relative to $(O/H)^{\text{tot}}$, following the prescriptions of Nicholls et al. (2017, their table 2). We update the values of $(N/H)^{\text{tot}}$ and $(C/H)^{\text{tot}}$ by randomly drawing from the chosen distributions of $(N/O)^{\text{tot}}$ and $(C/O)^{\text{tot}}$ (Section 2.2). In this approach, $(N/O)^{\text{tot}}$, $(C/O)^{\text{tot}}$, and $(O/H)^{\text{tot}}$ are therefore not correlated with one another, nor are they with the ionization parameter, which we draw fully independently. Any relation between

these quantities in galaxies should therefore arise from comparisons of models with observations. In fact, the results of Schaefer et al. (2020) show the importance of not imposing a fixed relation between, for example, $(N/O)^{\text{tot}}$ and $(O/H)^{\text{tot}}$ in emission-line studies of star-forming galaxies. We also perturb all other metal abundances around their initial values. Once a complete set of total abundances is obtained, we compute the associated metallicity $Z$, which we take to be the same for the ionizing star cluster.

We account for the depletion of all refractory elements on to dust grains according to the values in table 2 of Dopita et al. (2013), scaled in such a way as to match the randomly drawn dust-to-metal mass ratio, $\xi_d$. The resulting set of gas-phase abundances is used as input for CLOUDY. In the remainder of this paper, all abundances refer to the gas phase, except where indicated by a superscript 'tot'. Finally, the dust-to-gas ratio (DTG) entered in CLOUDY is determined as $\text{DTG/DTG}_{\odot} = \xi_d Z / (\xi_{d,\odot} Z_{\odot})$, with $\xi_{d,\odot} = 0.36$ and $Z_{\odot} = 0.01425$. For simplicity, we use the 'ism' dust type of CLOUDY with argument DTG/DTG$_{\odot}$, meaning that there is no strict coherence between element-by-element depletion and dust composition.

### 2.1.2 Stellar ionizing spectra

The models presented in this study rely on a set of ionizing spectral energy distributions computed using the Binary Population and Spectral Synthesis (BPASS) v2.2³ model of binary-star populations (Stanway & Eldridge 2018; Byrne et al. 2022). We adopt a fixed Chabrier (2003) initial mass function with lower and upper cut-offs 0.1 and 300 $M_{\odot}$, respectively.

We note that although some recent studies have suggested that X-ray-binary accretion discs may significantly contribute to the ionizing radiation in metal-poor star-forming galaxies (Schaerer, Fragos & Izotov 2019; Umeda et al. 2022; Katz et al. 2023), models calibrated on the observed average luminosity function of X-ray binaries in nearby galaxies do not support this hypothesis (Lecroq et al. 2024, see also Senchyna et al. 2020). We ignore this contribution in our analysis.

## 2.2 Production of H II-region models

The grid of H II-region models presented here is based on the set of adjustable parameters described in the previous sections. We draw values of the adjustable parameters from distributions defined in most cases by three variables: the minimum and the maximum values of a 'plateau' of constant probability, and the width $\sigma_w$ of the Gaussian function used to describe the drop of probability on both sides of this plateau. The values controlling the parameter distributions are listed in Table 1, where we indicate in each case the range and the width $\sigma_w$: some distributions have smooth edges, parametrized by the width $\sigma_w = 0.5$ of the Gaussian distribution defining the wings (Section 2.1). The value $\sigma_w = 0$ corresponds to a distribution strictly limited to the plateau range. We adopt quasi-log-uniform distributions for most parameters to effectively cover broad value ranges spanning multiple orders of magnitude, whereas the parameter HβFrac is sampled uniformly between 0.1 and 1.0 in steps of 0.1.

The $(C/O)^{\text{tot}}$ and $(N/O)^{\text{tot}}$ abundance ratios are allowed to vary in great proportions, while the other metals are only slightly perturbed (within ±0.1 dex) around the predefined values obtained from the Nicholls et al. (2017) prescription (Section 2.1.1). We allow this to explore how changes in the abundance of any single element can affect observables.

<footnote n="3">³https://bpass.auckland.ac.nz/</footnote>

<footnote n="23">²³MdB is a large database of photoionization models, with public access through the MySQL protocol since 2014 (http://3mdb.astro.unam.mx; see Morisset, Delgado-Inglada & Flores-Fajardo 2015).</footnote>
<table><caption>**Table 1.** Adjustable parameters used to generate the 3MdB/'HII_24' grid of CLOUDY H II region models, along with their distributions defined by the minimum and the maximum values of a 'plateau' of constant probability, and the width $\sigma_w$ of the Gaussian function used to describe the drop of probability on both sides of this plateau. A random shift $\Delta \log (X/H)^{\text{tot}}$ is applied to the total abundance of any metal X (other than C, N, and O) relative to the value defined in the adopted prescription (Section 2.1.1). The distributions refer to total abundances (i.e. including elements in both the gas and dust phases), the dust-to-metal mass ratio $\xi_d$ allowing transformation from these to the gas-phase values considered in all figures (see text for details).</caption><thead><tr><th>Parameter</th><th>min</th><th>max</th><th>σ<sub>w</sub></th></tr></thead><tbody><tr><td>log U</td><td>-4</td><td>-1</td><td>0.5</td></tr><tr><td>12 + log (O/H)<sup>tot</sup></td><td>6.5</td><td>9.5</td><td>0.5</td></tr><tr><td>log (N/O)<sup>tot</sup></td><td>-2</td><td>0</td><td>0.5</td></tr><tr><td>log (C/O)<sup>tot</sup></td><td>-1</td><td>0.5</td><td>0.5</td></tr><tr><td>Δ log (X/H)<sup>tot</sup></td><td>-0.1</td><td>0.1</td><td>0.0</td></tr><tr><td>t<sub>c</sub>/Myr</td><td>1</td><td>7</td><td>0.0</td></tr><tr><td>log f<sub>r</sub></td><td>-1.5</td><td>1.5</td><td>0.0</td></tr><tr><td>log(n<sub>H</sub>/cm<sup>-3</sup>)</td><td>1</td><td>4</td><td>0.5</td></tr><tr><td>log ε</td><td>-2</td><td>0</td><td>0</td></tr><tr><td>ξ<sub>d</sub></td><td>0.15</td><td>0.55</td><td>0.0</td></tr><tr><td>Hβ Frac</td><td colspan="3">From 0.1 to 1.0 in steps of 0.1</td></tr></tbody></table>

The 3MdB data base currently holds, under the reference 'HII_24', about 3.5 million models drawn from the distributions listed in Table 1. These were obtained from 400 000 runs of CLOUDY, as each run produces up to 10 entries in the data base once the optical depth cuts are applied (Section 2.1). The models include both families of BPS and BPB ionizing spectra, corresponding to BPASS single- and binary-star populations respectively (Section 2.1.2). We primarily focus in this paper on results obtained with 2001 534 BPB models.

In Fig. 1, we plot the distributions of various quantities obtained for 1 million BPB models with $t_c < 7$ Myr and $H\beta$ Frac > 0.65 (corresponding to $f_{esc} \lesssim 0.35$; Section 2.1). As expected, distributions with $\sigma_w = 0$ (e.g. for the age $t_c$) exhibit sharp edges, and those with $\sigma_w = 0.5$ (e.g. for $\log n_H$) exhibit smooth edges. The distributions of element abundances are in general not entirely regular, as they are defined indirectly. Indeed, the distributions in Table 1 pertain to total (gas + dust phase) abundances, while those shown in Fig. 1 refer to gas-phase abundances, affected by the dust-to-metal mass ratio $\xi_d$ drawn on a case-by-case basis.

The small drop between -5.5 and -4.5 in the log(S/H) distribution results from the change in slope of the relation linking (S/H)tot to (O/H)tot according to Nicholls et al. (2017, see their fig. 8). There is no similar drop in the distributions of Ar/H and Ne/H, as the relations defining the abundances of these elements relative to oxygen do not exhibit any change in slope (table 2 of Nicholls et al. 2017). We note the small tails of $H\beta$ Frac-distribution spikes below the discrete values drawn from Table 1. This is because we do not interpolate between 'zones' (radial steps) of computation defined by CLOUDY and only integrate the $H\beta$ luminosity out to the zone ending closest to the selected $H\beta$ Frac. Similarly, the log U distribution shows a small tail towards large values, since the distribution in Table 1 pertains to the radiation-bounded model, while the volume-averaged ionization parameter tends to increase towards the interior of the nebula and is then higher in the matter-bounded models ($H\beta$ Frac < 1.0).

The simplest way to populate a multidimensional grid is to assign parameters in each dimension to regularly spaced values, in either linear or logarithmic scales, as used by e.g. Morisset et al. (2016), Vale Asari et al. (2016), Galliano (2018), Ramambason et al. (2022), and Smirnov-Pinchukov et al. (2022). For extremely large

grids, a more practical approach involves using randomly distributed independent parameters, which offers an efficient way to populate data sets to train ML algorithms (e.g. Bron et al. 2021). In future iterations of our method to generate a grid of H II regions, we plan to adopt more efficient distributions of input parameters. This will help mitigate random clustering of data points in certain regions of the parameter space, which could introduce biases into the training of the ANN. Specifically, we aim to explore methods such as Latin Hypercube Sampling and quasi-random sampling techniques such as Sobol and Halton sequences, as used by, e.g. Mootoovaloo et al. (2022), Bird et al. (2019), and Reza et al. (2022), as well as the method described by Korobov (1959) and applied by, e.g. DeRose et al. (2022) and Ramirez et al. (2024). With such approaches, the plateau-like parameter distributions with smoothly decreasing edges used in this work will no longer be preserved.

## 3 REGRESSION BASED ON MACHINE LEARNING TO INTERPOLATE NEBULAR-EMISSION MODELS

In this section, we use ML to interpolate unlimited amounts of models from the grid presented in Section 2. We present the regression model, evaluate its performance, and explore a few applications by examining the influence of the main adjustable parameters on selected observables. Our approach consists in creating, training, and testing a regression ANN, here managed through the AI4NEB library<fnref n="4" /> to predict the properties of models with any given input parameters. The neural network counts three hidden layers with 256, 512, and 256 activation cells, respectively. All activation functions are of the Gaussian Error Linear Unit (GELU) type (Hendrycks & Gimpel 2016). A detailed description of the ANN is provided in Appendix A.

### 3.1 Regression model

We use 800 000 models for the training set (a random selection of 80 per cent of the models in Fig. 1) to train an ANN receiving 13 parameters as input,

(i) log U

(ii) log $f_r$

(iii) log $t_c$

(iv) log(O/H)

(v) log(C/O)

(vi) log(N/O)

(vii) log(Ne/O)

(viii) log(S/O)

(ix) log(Ar/O)

(x) log $n_H$

(xi) $H\beta$ Frac

(xii) $\xi_d$

(xiii) $\epsilon$,

and predict, as output, 14 line-luminosity ratios, the electronic temperature of the ionized gas, and the $H\alpha$ luminosity (all in logarithmic scale),

(i) $H\alpha/H\beta$

(ii) [O III] $\lambda$5007/H$\beta$

<footnote n="4">⁴AI4NEB allows the trained network to be saved, along with the scaling function (used to normalize the input data and the training and test sets), as described by Morisset et al. (in preparation). See https://github.com/Morisset/AI4neb.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Distributions of selected input parameters of 1 million CLOUDY models with $t_c < 7$ Myr and $H\beta$Frac > 0.65 produced using the distributions in Table 1. Also shown in the bottom right-hand corner is the distribution of output Hα luminosities.</figcaption></figure>

(iii) [N II] λ6584/Hβ

(iv) [O II] λ3726, 3729/Hβ (hereafter [O II] λ3727/Hβ)

(v) [O III] λ4363/λ5007

(vi) [N II] λ5755/λ6584

(vii) [S II] λ6716, 6731/Hβ (hereafter [S II] λ6725/Hβ)

(viii) [S III] λ9530/Hβ

(ix) [Ar III] λ7135/Hβ

(x) [Ne III] λ3868/Hβ

(xi) [S II] λ6716/λ6731

(xii) He I λ5876/Hβ

(xiii) [S III] λ6312/λ9530

(xiv) [O II] λ7320, 7330/Hβ (hereafter [O II] λ7325/Hβ)

(xv) $T_c$

(xvi) $H\alpha$.

The list of inputs is in practice redundant, as, for example, S/O is fully determined once O/H and $\xi_d$ are known. Such redundancy is commonly employed to enhance the performance of regression models by providing extra information strongly related with the desired outputs – in this case, the luminosities of S emission lines (see e.g. Ho 2019, where redundant line ratios and sums are used).

We have included the electronic temperature as an output variable, because it is treated as such (and not as an input parameter) in our CLOUDY models and will enable exploration of the influence of temperature fluctuations on the integrated line emission from a collection of H II regions (Section 4). Also, the Hα luminosity is a key output quantity to identify models that may result from unlikely combinations of randomly drawn input parameters.

## 3.2 Performance of the regression model

Fig. 2 shows the distributions of the differences between the true values of all 14 line ratios and the electronic temperature (all in logarithmic scale) and the values of these quantities retrieved from the regression model. To evaluate the quality of the ANN, we focus on the subset of 200 000 test-set points projecting into a specific triangular region of the [O III] λ5007/Hβ-versus-[N II] λ6584/Hα 'BPT' diagram (Baldwin, Phillips & Terlevich 1981). This region corresponds to the location of the models used in subsequent steps and is defined by the following three conditions: $[O \text{ III}] \lambda 5007/H\beta < 10$; $[N \text{ II}] \lambda 6584/H\alpha < 1$; and $\log([O \text{ III}] \lambda 5007/H\beta) > -1.1 \log([N \text{ II}] \lambda 6584/H\alpha) - 2.6$. Applying this filter selects 61 110 data points from the 200 000 available in the test set. Distributions are mapped as a function of oxygen abundance, with orange lines indicating the 5th–95th percentile range. We also indicate, at the top of each panel, the global standard deviation of the difference between true and retrieved values. To complement Fig. 2, we show in Fig. 3 the difference between true and predicted $\log([O \text{ III}] \lambda 5007/H\beta)$ ratio as a function of each of the 13 input parameters of the regression model. Finally, in Fig. 4, we compare true (stars) and predicted (circles) line ratios for a random subset of 40 models in the BPT diagram, colour-coded according to oxygen abundance (left) and ionization parameter (right).

Figs 2–4 together demonstrate the high quality of the regression model. The standard deviation in Fig. 2 is less than 0.02 dex for all quantities, in stark contrast to the scatter in corresponding observed values, which can exceed 1 dex in data sets such as the SDSS
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Density maps of the difference between true and ANN-predicted values (both in logarithm) of 14 emission-line ratios, the electronic temperature, and the absolute Hα luminosity, plotted against 12 + log (O/H), for the 20 per cent of models from Fig. 1 not included in the ANN training set (darker gray corresponding to higher model density). In each panel, the orange lines show the 5th–95th percentile range of the distribution as a function of 12 + log (O/H), while the standard deviation of the differences between true and predicted logarithmic values is indicated at the top.</figcaption></figure>

(Sloan Digital Sky Survey, see Fig. 5 discussed in the next section). Moreover, the errors in the predictions are relatively symmetric, being equally distributed between positive or negative differences. This reinforces our confidence in the regression model, as small symmetric errors can compensate for each other when exploring trends with a large number of predicted models. In Fig. 4, true and predicted models are also very close to each other, with no systematic effect. Fig. 2 shows a clear trend of increasing errors towards large O/H: at metallicities above 12 + log (O/H) ≈ 9, the regression model starts to lose accuracy. At such high metallicities, the low electron temperature amplifies the impact of small parameter variations, which can be more difficult to track with the ANN model. Such metallicities are barely reached in the applications presented in the remainder of this paper.

We conclude that the regression ANN presented in this section provides an efficient means to predict the observable properties of individual H II regions. It allows us to generate a large number of photoionization models and explore the influence of input parameters on emission-line predictions much faster than would be required with CLOUDY.

Our ANN can reliably predict (Figs 2–4 and Table A1) the electronic temperature, Hα luminosity, and strengths of 14 line ratios for 10⁶ combinations of the 13 input parameters, corresponding to as many H II regions, in less than 0.2 s using a single Graphics Processing Unit (GPU), i.e. nine orders of magnitude faster than

running CLOUDY on a single Central Processing Unit (CPU). Details on the performance of the model predictions, including memory and time usage for different grid sizes, are provided in Table 2, comparing CPU and GPU configurations.

We note that the ANN can also be used as part of an algorithm to search for the parameters that best reproduce observations of a given H II region [such as Markov chain Monte Carlo (MCMC) and genetic algorithms]. This will be the subject of a forthcoming paper.

## 3.3 Exploring the influence of input parameters on emission-line predictions for H II regions

With the ANN developed to enable fast modelling of H II regions, we can efficiently explore how basic relations between input parameters, such as log *U* versus O/H, impact the predicted observables in classical diagnostic diagrams. This analysis has didactic value, as the results can be directly compared with survey data, to provide insight into the validity and quality of particular assumptions. In this section, we explore some of these relations.

The four central panels of Fig. 5 show the locations of 10 000 models with both N/O and *U* tied to O/H (according to the relations drawn in the leftmost panels) in four line-ratio diagrams: the [O III] λ5007/Hβ-versus-[N II] λ6584/Hα (BPT) diagram; the [O III] λ5007/Hβ-versus-[S II] λ6725/Hα and [O III] λ5007/Hβ-versus-[O II] λ3727/Hα diagrams (Veilleux & Osterbrock 1987,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Analogue of Fig. 2 for the logarithmic difference between true and predicted [O III] λ5007/Hβ ratio as a function of each of the 13 input parameters of the ANN of Section 3.2.</figcaption></figure>

hereafter VO87); and the $P_2$-versus-$P_1$ diagram defined by Ji & Yan (2022) by combining $H\alpha$, [N II], [S II], and [O III] lines.<fnref n="5" /> The models have fixed ionizing-cluster age, $t_c = 3$ Myr, hydrogen density, $n_H = 30$ cm⁻³, form factor, $f_r = 1.0$, HβFrac = 0.95, and dust-to-metal mass ratio $\xi_d = 0.35$. Also shown in the four central panels of Fig. 5 are the observed properties of 28 311 star-forming (i.e. BPTCLASS = 1) SDSS galaxies from the DR8 MPA-JHU catalogue,<fnref n="6" /> selected to have signal-to-noise ratio S/N > 10 in all strong lines of interest and S/N > 3 in He I λ5876 and [Ne III] λ3869 (orange background). In the rightmost panels, we compare the models to observations (in black) of 603 H II regions compiled by Marino et al. (2013) in two diagrams showing the temperature-sensitive [O III] λ4363/5007 ratio against 12 + log (O/H) (top) and the ([O II] λ3727 + [O III] λ5007)/Hβ ratio (bottom).

In all panels of Fig. 5, the models follow a narrow sequence set by the unique relation between O/H, N/O, and U chosen for the purpose of this illustration. The sequence overlaps with at least some observations in all line-ratio diagrams, but does not reflect the observed scatter. In the $P_1-P_2$ diagram, unlike in other diagrams, the models do not reproduce at all the observational sequence, particularly along the $P_2$ axis. The sequence is captured when N/O or log U are allowed to deviate from the fixed relation, as discussed below.

We explore the effect of allowing for variations in log(N/O) and log U at fixed 12 + log (O/H) in Fig. 6. The top eight panels show the results for 10 000 models as in Fig. 5, but drawing N/O uniformly withing ±0.3 dex of the relation in the top-left panel of Fig. 5.

The bottom eight panels of Fig. 6 show the results for 10 000 models as in Fig. 5, but drawing U uniformly withing ±0.5 dex of the relation in the bottom-left panel of Fig. 5. For completeness, in these panels we also draw 1000 models for which the ionizing parameter is set to increase with oxygen abundance, as found by Ji & Yan (2022, their fig. 10) for extragalactic H II regions selected from the SDSS-MANGA survey (Mapping Nearby Galaxies at Apache Point Observatory), which we parametrize roughly as $log U = -7.2 + 0.5[12 + log (O/H)]$ (in red).

Fig. 6 illustrates how modest changes in the adopted dependence of N/O and U on O/H can significantly affect the predicted line ratios. The influence of changes in N/O at fixed O/H is strong on ratios involving [N II] lines, but negligible on other line ratios. In contrast, changes in the dependence of U on O/H significantly affect all line ratios, including the $T_c$-sensitive [O III] λ4363/λ5007 ratio. Such variations allow the models to cover a larger portion of the observational sequences in the various diagrams. We note that the model where log U increases with 12 + log (O/H) leads to predictions far from the observations at low metallicity. This positive correlation identified by Ji & Yan (2022) may result from their use of a fixed relation between N/O and O/H, following the prescription of Dopita et al. (2013).

In Appendix B, we explore the effect on line ratios of changing the other adjustable parameters of the models, namely ionizing-cluster age, hydrogen density, HβFrac, and dust-to-metal mass ratio.

## 4 BUILDING THE NEBULAR EMISSION FROM COMPOSITE STAR-FORMING GALAXIES

The ANN presented in the previous section allows us to calculate in record time the emission-line properties of large numbers of H II regions with arbitrary physical parameters. In this section, we exploit this tool to build synthetic star-forming galaxies by efficiently

<footnote n="5">⁵$P_1 = 0.63[\text{N II}] \lambda6584/\text{H}\alpha + 0.51[\text{S II}] \lambda6725/\text{H}\alpha + 0.59[\text{O III}] \lambda5007/\text{H}\beta$ and $P_2 = -0.63[\text{N II}] \lambda6584/\text{H}\alpha + 0.78[\text{S II}] \lambda6725/\text{H}\alpha$.</footnote>

<footnote n="6">⁶Max Planck Institute for Astrophysics and the Johns Hopkins University
Techniques: https://www.sdss4.org/dr17/spectro/galaxy mpajhu/</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** True (stars) and predicted (circles) line ratios in the [O III] λ5007/Hβ-versus-[N II] λ6584/Hα (BPT) diagram, for a random subset of 40 models from the sample shown in Fig. 2. The models are colour-coded according to oxygen abundance (left) and ionization parameter (right).</figcaption></figure>

generating and assembling large populations of H II regions with different parameters, including individual gas-covering factors. These regions are individually attenuated by dust, and a component of diffuse ionized gas (DIG) is added to the composite galaxy. Then, we apply a global correction for dust attenuation to the emission-line luminosities of these synthetic galaxies using the Hα/Hβ ratio, following the approach typically employed for observed galaxies. We describe the detailed steps of this process in the following sections.

## 4.1 Distribution of H II-region parameters

We compute the nebular emission from composite star-forming galaxies by summing the emission-line luminosities from a collection of H II regions. The emission lines considered here are all optically thin and are affected only by dust attenuation, which is accounted for by assigning a randomly determined A$_{V}$ value to each H II region (Section 4.2).

For each galaxy model, we start by drawing a ‘baseline’ galactic oxygen abundance, (O/H)gal, around which we draw the oxygen abundances of individual H II regions, O/H. We write

$$12 + \log (\text{O/H})_{\text{gal}} = 8 + \varrho(\delta_{\text{O/H}}), \quad (1)$$

where $\varrho(\delta_{\text{O/H}})$ is a uniformly distributed random function taking values between $-\delta_{\text{O/H}}$ and $+\delta_{\text{O/H}}$, and

$$\log(\text{O/H}) = \log(\text{O/H})_{\text{gal}} + \varrho(\sigma_{\text{O/H}}). \quad (2)$$

Hence, the quantity $\delta_{\text{O/H}}$ defines the departure of the baseline O abundance from the canonical $12 + \log (\text{O/H})_{\text{gal}} = 8$ value, while the quantity $\sigma_{\text{O/H}}$ defines the random variation of O/H around this baseline abundance from one H II region to another, in a given galaxy.

We define other nebular parameters in a similar way, by drawing baseline galactic values through standard relations of the parameter of interest to O/H, and H II-region values by sampling around these baseline galactic values. Specifically, we assume that the ionization parameter scales roughly as the inverse of metallicity (Carton et al. 2017), adopting a canonical value $\log U_{\text{gal}} = -2.5$ at the canonical O abundance $12 + \log (\text{O/H})_{\text{gal}} = 8$ (equation 1). We thus write

$$\log U = -2.5 - [\log(\text{O/H}) + 4] + \varrho(\delta_{\text{U}}) + \varrho(\sigma_{\text{U}}). \quad (3)$$

For the dependence of C/O, N/O, Ne/O, S/O, and Ar/O on O/H, we start from the baseline relations derived by Nicholls et al. (2017, their table 2) for total abundances, which we convert into relations for gas-phase abundances via the dust-to-metal mass ratio $\xi_d$ tied to the depletion factors of Dopita et al. (2013, see Section 2.1.1 above). We achieve this by computing, for a grid of (O/H)tot and $\xi_d$ sampling the full ranges of these quantities in Table 1, the corresponding gas-phase abundances O/H, C/O, N/O, Ne/O, S/O, and Ar/O. We use this data set to train a small ANN (called abund-ANN) to interpolate C/O, N/O, Ne/O, S/O, and Ar/O based on O/H and $\xi_d$.<fnref n="7" /> With this framework, for each galaxy for which a baseline galactic O/H was selected from equation (2), we select a baseline dust-to-metal mass ratio uniformly in the range $0.2 \le \xi_d \le 0.5$ and use the abund-ANN to determine the corresponding abundance ratios (X/O)fit, where X stands for N, Ne, S, or Ar. By analogy with equations (1) and (2), we define baseline abundance ratios for this galaxy as

$$\log (\text{X/O})_{\text{gal}} = \log (\text{X/O})_{\text{fit}} + \varrho(\delta_{\text{X/O}}), \quad (4)$$

and the abundance ratios of individual H II regions as

$$\log(\text{X/O}) = \log(\text{X/O})_{\text{gal}} + \varrho(\sigma_{\text{X/O}}). \quad (5)$$

We note that the quantities $\varrho(\delta_{\text{U}})$ and $\varrho(\delta_{\text{N/O}})$ in equations (3) and (4) above provide a means of investigating small changes in the adopted dependence of N/O and U on O/H, whose importance was noted in Section 3.3 (Fig. 6). The quantity $\varrho(\delta_{\text{N/O}})$ can be assimilated to the nitrogen excess factor considered by Schaefer et al. (2020) and denoted by $\xi$ in their nomenclature.

In Table 3, we list the spread in baseline galactic values (second column), and that in individual H II-region values around a given baseline value (third column), for the oxygen abundance, ionization parameter, and the different metal-abundance ratios explored here (top seven lines). To probe the influence of star formation history, we randomly define for each galaxy the minimum and maximum H II-region ages, $t_c^{\text{min}}$ and $t_c^{\text{max}}$, within the range of cluster ages used to train the ANN (Table 1). We then randomly draw the ages of all the galaxy’s H II regions between $t_c^{\text{min}}$ and $t_c^{\text{max}}$. This allows us to include galaxies containing only very young bursts, only older bursts, or a broad distribution of burst ages, as might be the case, for example, with constant star formation. For the form factor $f_r$, hydrogen density

<footnote n="7">⁷The abund-ANN is developed using the Keras library (https://github.com/keras-team/keras) running on the TENSORFLOW back end (http://tensorflow.org/). It consists of two hidden layers, each with 50 neurons using GELU activations (Hendrycks & Gimpel 2016), trained with the Adam optimizer (Kingma & Ba 2014) for 200 epochs on a training set of 32 000 models. The input is of dimension two and the output is of dimension five. The standard deviations of the differences between the predictions and the test values are of the order of 0.002. The abund-ANN is managed through the AI4NEB library.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Properties of 10 000 regression-ANN models with both N/O and U tied to O/H (according to the relations shown in the leftmost panels) in six diagrams: the [O III] λ5007/Hβ-versus-[N II] λ6584/Hα (BPT) diagram; the [O III] λ5007/Hβ-versus-[S II] λ6725/Hα and [O III] λ5007/Hβ-versus-[O II] λ3727/Hα diagrams (VO87); the P₂-versus-P₁ diagram defined by Ji & Yan (2022); and two diagrams showing the temperature-sensitive [O III] λ4363/5007 ratio against 12 + log (O/H) (top right) and the ([O II] λ3727 + [O III] λ5007)/Hβ ratio (bottom right). The models have fixed ionizing-cluster age, t_c = 3 Myr, hydrogen density, n_H = 30 cm⁻³, form factor, f_t = 1.0, HβFrac = 0.95, and dust-to-metal mass ratio ξ_d = 0.35. The models are colour-coded according to oxygen abundance (see scale on the right). The orange symbols (in the four central panels) show observations of SDSS star-forming galaxies, and the black symbols (in the rightmost panels) show observations of H II regions selected from the SDSS-MANGA survey (see text for more details).</figcaption></figure>

<table><caption>**Table 2.** Performance of the regression ANN in predicting emission-line properties of H II regions across various CPU and GPU configurations. The first column lists the number of models run by the ANN, while the second and third columns indicate the number of CPUs and GPUs used. The fourth column shows the peak Random-Access Memory (RAM) usage for CPUs, no data being given when GPU is used. The last column provides execution time.</caption><thead><tr><th>N models</th><th>CPU<sup>a</sup></th><th>GPU<sup>b</sup></th><th>Max memory used (GB)</th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="6">2.0 × 10<sup>5</sup></td><td>1</td><td>0</td><td>2.8</td><td>2.08</td></tr><tr><td>12</td><td>0</td><td>2.8</td><td>0.27</td></tr><tr><td>12<sup>c</sup></td><td>0</td><td>2.8</td><td>0.65</td></tr><tr><td>12<sup>d</sup></td><td>0</td><td>2.8</td><td>1.8</td></tr><tr><td>0</td><td>1</td><td>—</td><td>0.036</td></tr><tr><td>1</td><td>0</td><td>16</td><td>12</td></tr><tr><td rowspan="6">1.6 × 10<sup>6</sup></td><td>12</td><td>0</td><td>16</td><td>1.75</td></tr><tr><td>12<sup>c</sup></td><td>0</td><td>16</td><td>9.6</td></tr><tr><td>12<sup>d</sup></td><td>0</td><td>16</td><td>36<sup>e</sup></td></tr><tr><td>60</td><td>0</td><td>16</td><td>1.15</td></tr><tr><td>0</td><td>1</td><td>—</td><td>0.26</td></tr><tr><td>1</td><td>0</td><td>123</td><td>131</td></tr><tr><td rowspan="3">1.3 × 10<sup>7</sup></td><td>12</td><td>0</td><td>123</td><td>13.2</td></tr><tr><td>30</td><td>0</td><td>123</td><td>8.8</td></tr><tr><td>60</td><td>0</td><td>123</td><td>8.8</td></tr></tbody></table>

$n_H$, volume-filling factor of the gas $\epsilon$, and dust-to-metal mass ratio $\xi_d$, we sample over the full ranges used to train the ANN, while we restrict HβFrac to values in the range 0.7–1.0, corresponding to (relatively) radiation-bounded H II regions.

To simulate star-forming galaxies as realistically as possible, we must also specify a luminosity distribution of H II regions such that the observed Hα luminosity function is reproduced. First, we account for the fact that, in addition to being potentially matter-bounded, H II regions can also lose ionizing photons through holes resulting from incomplete covering of the ionizing source by the surrounding gas. Following Heckman et al. (2011), we parametrize such a ‘picket-fence’ model in terms of the fraction $f_c = \Omega/4\pi$ of solid angle over which the gas covers the ionizing star cluster, which we randomly draw between 0.1 and 1.0 for each H II region. The effect is that all line luminosities are reduced by a factor $f_c$, without affecting the line ratios.

Then, for the Hα luminosity function of H II regions, we adopt a power-law distribution LF(LHα) ∝ LHα$\alpha_{LF}$, where LF(LHα)/dLHα is the number of regions with luminosities between LHα and LHα + dLHα. We fix the slope at $\alpha_{LF} = -1.9$, intermediate between the values −2.12 ± 0.03 and −1.73 ± 0.15 found from spatially resolved spectroscopy of ~4300 H II regions in the nearby, face-on spiral galaxy NGC 628 (Rousseau-Nepton et al. 2018)<fnref n="8" /> and ~23 000 H II regions in a sample of 19 nearby, low-inclination spiral galaxies (Santoro et al. 2022).

In practice, we build the emission-line properties of a composite star-forming galaxy by (i) drawing the physical parameters of a large sample of H II regions using the procedure described above (Table 3); (ii) using the ANN to compute the associated emission-line properties of these H II regions (Section 3.1); (iii) reducing the emission-line luminosities by the gas-covering factor selected for each H II region; and (iv) including individual H II regions in the galaxy emission-line budget or rejecting them on the basis of their Hα luminosity,

<footnote><sup>a</sup>AMD EPYC 7313 processor, unless otherwise specified.</footnote>

<footnote><sup>b</sup>NVIDIA RTX A5000 graphics card.</footnote>

<footnote><sup>c</sup>M2 Pro processor on a MacMini.</footnote>

<footnote><sup>d</sup>Intel Core i7 2.6 GHz processor on a MacBook Pro.</footnote>

<footnote><sup>e</sup>In this configuration, the RAM was filled, triggering a swap process, which increased the execution time.</footnote>

<footnote n="8">⁸The numerical values of the quantity $\alpha$ quoted in Rousseau-Nepton et al. (2018) refer to the slope of the logarithmic Hα luminosity function, $\alpha = \alpha_{LF} + 1$ (Rousseau-Nepton, private communication).</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 6. The top eight panels show the analogue of Fig. 5 for the same 10 000 models, but drawing N/O uniformly withing ±0.3 dex of the relation in the top-left panel of Fig. 5 (as indicated by the colour scale on the right). The bottom eight panels show the analogue of Fig. 5 for the same 10 000 models, but drawing U uniformly withing ±0.5 dex of the relation in the bottom-left panel of Fig. 5 (as indicated by the colour scale on the right). Also shown in red are the results for 10 000 models with log U = −7.2 + 0.5[12 + log (O/H)] (following Ji & Yan 2022). The data are the same as in Fig. 5.</figcaption></figure>

as determined by the target luminosity function. We accomplish this
last step by associating with each H II region a chance

$$
C = \begin{cases}
\frac{(\mathrm{LH}\alpha / \mathrm{LH}\alpha_{\mathrm{min}})^{\alpha_{\mathrm{LF}}}}{\nu(0,1)} & \text{if } \mathrm{LH}\alpha_{\mathrm{min}} < \mathrm{LH}\alpha < \mathrm{LH}\alpha_{\mathrm{max}}, \quad (6) \\
0 & \text{otherwise,}
\end{cases}
$$

where $\nu(0, 1)$ is a random variable drawn for each H II region from
a uniform distribution over the interval from 0 to 1. The exact values
of $\mathrm{LH}\alpha_{\mathrm{min}}$ and $\mathrm{LH}\alpha_{\mathrm{max}}$ have little influence on our results. We
adopt $\mathrm{LH}\alpha_{\mathrm{min}} = 1 \times 10^{36} \mathrm{erg}\,\mathrm{s}^{-1}$ and $\mathrm{LH}\alpha_{\mathrm{max}} = 3 \times 10^{40} \mathrm{erg}\,\mathrm{s}^{-1}$
(Rousseau-Nepton et al. 2018; Santoro et al. 2022). An H II region is
included in the galaxy emission-line budget only if the chance given
by equation (6) exceeds unity, i.e. $C > 1$.

We compute in this way the emission-line properties of 250 000
composite star-forming galaxies containing at least N_HII,min = 100
individual H II regions. Each galaxy is built by drawing an initial
number of 50 000 H II regions, the majority of which are rejected
based on the above criterion on C, so that model galaxies in
practice contain up to about 3000 H II regions. The low efficiency
of the process is compensated for by the speed of the ANN. A
higher efficiency could be achieved by choosing a different set
of input parameters for the regression model (Section 3.1), for
example, by replacing the ionization parameter and form factor
with the ionizing-photon rate and inner radius of the H II region.
The disadvantage is that this would not allow exploration of the
influence of the relation between ionization parameter and oxygen
<table><caption>**Table 3.** Adjustable parameters used to generate composite star-forming galaxies (see Section 4.1 for details).</caption><thead><tr><th>Parameter</th><th>Galaxy spread</th><th>H II-region spread</th><th>Description</th></tr></thead><tbody><tr><td>log (O/H)</td><td>δ<sub>O/H</sub> = 1.2</td><td>σ<sub>O/H</sub> = 0.15</td><td>Oxygen abundance (see equations 1 and 2 for canonical value)</td></tr><tr><td>log U</td><td>δ<sub>U</sub> = 0.75</td><td>σ<sub>U</sub> = 0.5</td><td>Ionization parameter (canonical value from equation 3)</td></tr><tr><td>ξ<sub>d</sub></td><td>[0.20, 0.50]</td><td>σ<sub>ξd</sub> = 0.05</td><td>Dust-to-metal mass ratio</td></tr><tr><td>log (C/O)</td><td>δ<sub>C/O</sub> = 0.3</td><td>σ<sub>C/O</sub> = 0.2</td><td>C/O abundance ratio (canonical value tied to O/H and ξ<sub>d</sub>; Section 4.1)</td></tr><tr><td>log (N/O)</td><td>δ<sub>N/O</sub> = 0.3</td><td>σ<sub>N/O</sub> = 0.2</td><td>N/O abundance ratio (canonical value tied to O/H and ξ<sub>d</sub>; Section 4.1)</td></tr><tr><td>log (S/O)</td><td>δ<sub>S/O</sub> = 0.3</td><td>σ<sub>S/O</sub> = 0.2</td><td>S/O abundance ratio (canonical value tied to O/H and ξ<sub>d</sub>; Section 4.1)</td></tr><tr><td>log (Ne/O)</td><td>δ<sub>Ne/O</sub> = 0.3</td><td>σ<sub>Ne/O</sub> = 0.2</td><td>N/O abundance ratio (canonical value tied to O/H and ξ<sub>d</sub>; Section 4.1)</td></tr><tr><td>log (Ar/O)</td><td>δ<sub>Ar/O</sub> = 0.3</td><td>σ<sub>Ar/O</sub> = 0.2</td><td>Ar/O abundance ratio (canonical value tied to O/H and ξ<sub>d</sub>; Section 4.1)</td></tr><tr><td>Hβ/Frac</td><td>[0.75, 0.95]</td><td>σ<sub>HbFrac</sub> = 0.05</td><td>Hβ-luminosity ratio of truncated versus radiation-bounded models</td></tr><tr><td>DIGFrac</td><td>[0.0, 0.4]</td><td>N/A</td><td>Hβ-luminosity ratio of DIG versus H II regions</td></tr><tr><td>t<sup>min</sup><sub>c</sub>/Myr</td><td>[1, 7]</td><td>N/A</td><td>Minimum ionizing star-cluster age</td></tr><tr><td>t<sup>max</sup><sub>c</sub>/Myr</td><td>[t<sup>min</sup><sub>c</sub>, 7]</td><td>N/A</td><td>Maximum ionizing star-cluster age</td></tr><tr><td>t<sub>c</sub>/Myr</td><td>N/A</td><td>[t<sup>min</sup><sub>c</sub>, t<sup>max</sup><sub>c</sub>]</td><td>H II-region age</td></tr><tr><td>log(n<sub>H</sub>/cm<sup>-3</sup>)</td><td>N/A</td><td>[0.5, 3]</td><td>Hydrogen density</td></tr><tr><td>log ε</td><td>N/A</td><td>[-2, 0]</td><td>Volume-filling factor of the gas</td></tr><tr><td>A<sub>V</sub>/mag</td><td>N/A</td><td>[0, 2]</td><td>V-band dust attenuation</td></tr><tr><td>log f<sub>r</sub></td><td>N/A</td><td>[-1.5, 1.5]</td><td>Geometric form factor</td></tr><tr><td>f<sub>c</sub></td><td>N/A</td><td>[0.1, 1.0]</td><td>Gas-covering factor</td></tr><tr><td>α<sub>LF</sub></td><td>N/A</td><td>N/A</td><td>Slope of the Hα luminosity function of H II regions, fixed at -1.9</td></tr><tr><td>N<sub>HII,min</sub></td><td>N/A</td><td>N/A</td><td>Minimum number of H II regions, fixed at 100</td></tr></tbody></table>

abundance (equation 3) on the diagnostic diagrams described in Section 3.3.

## 4.2 Components linked to molecular clouds and the diffuse ISM

We combine the emission from H II regions in composite star-forming galaxies with components linked to molecular clouds and the diffuse ISM.

To reflect the diversity of molecular-cloud environments in which stars form, we attenuate the emission from each individual H II region using the Galactic extinction curve of Fitzpatrick (1999) with $R_V = 3.1$, scaled by a V-band dust attenuation, $A_V$, drawn randomly between 0 and 2 (Table 3). This is justified by the fact that, while $A_V$ for a given region can be influenced by its size and metallicity, it also strongly depends on the presence of dust clouds along the line of sight, which is further affected by the orientation of the host galaxy and the position of the region within it. The specific wavelength dependence of dust attenuation has a negligible influence on the results presented in this study.

We also include a component of DIG, which is an important ingredient of nebular emission in star-forming galaxies (e.g. Sanders et al. 2017; Stasinska 2019). For this, we call on the CLOUDY calculations of Martínez-Paredes et al. (2023), which describe the emission of interstellar gas photoionized by populations of HOLMES (see Flores-Fajardo et al. 2011). These stars are expected to be the main sources of ionizing photons producing the DIG (Rautio et al. 2022; Postnikova & Bizyaev 2023). We do not consider here the potential (presumably minor) contribution from leakage of ionizing photons from H II regions (e.g. McClymont et al. 2024).

We use a regression ANN, called DIG-ANN, similar to that described in Section 3 above, to interpolate the DIG emission-line luminosities for each composite galaxy from the Martínez-Paredes et al. (2023) model grid. We add a fraction DIGFrac of the emission-line luminosities obtained in this way, drawn uniformly in the range $0.0 \le DIGFrac \le 0.4$, to the attenuated line emission from H II regions, to obtain the total line emission of the composite galaxy. The DIG-ANN and its training set are described in detail in Appendix C.

Then, we correct the total line emission of the composite galaxy based on the departure of the integrated $H\alpha/H\beta$ ratio from the

reference dust-free value of 2.86 (appropriate for a gas of density $100 \text{ cm}^{-3}$ and temperature $10\,000 \text{ K}$; Osterbrock & Ferland 2006), using the Galactic extinction curve of Fitzpatrick (1999), as would be done for an observed galaxy. We refer to the integrated, dust-corrected $H\alpha$ luminosity of the galaxy computed in this way as $LH\alpha_{gal}$.

As an example, we show in Fig. 7 the distributions of individual parameters of the 1616 H II regions constituting one such composite star-forming galaxy, for which the baseline ionization parameter and abundances (Section 4.1) are indicated by vertical dashed lines.

## 5 RESULTS

In this section, we explore the integrated properties of the 250 000 composite star-forming galaxies built as described in the previous section. We investigate the location of these galaxies in standard line-ratio diagnostic diagrams. Then, to illustrate the usefulness of this sample to interpret the integrated emission-line properties of composite star-forming galaxies, we examine the extent to which determinations of the oxygen abundance based on the integrated line emission of a galaxy differ from the average abundance of its constituent H II regions.

### 5.1 Integrated properties of composite star-forming galaxies

We show in Fig. 8 the distributions of integrated properties of the 250 000 composite star-forming galaxies built as described in the previous section. These include in particular the number of constituent H II regions ($N_{HII}$) and the integrated, dust-corrected $H\alpha$ luminosity ($LH\alpha_{gal}$), which incorporates the DIG contribution (Section 4.2). Also shown are average H II-region parameters, computed by weighting the parameters of individual H II regions by the intrinsic $H\alpha$ luminosities of these regions (LH$\alpha$; Section 4). In what follows, we refer to these average parameters as the ‘$H\alpha$-weighted H II-region parameters’.

The distributions of the various quantities in Fig. 8 can be understood from the way in which the parameters are sampled in Table 3. For example, the relation between the baseline $12 + \log (O/H)$ and $\log U$ (equation 3), together with the upper limit on LH$\alpha$ (equation 6)
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Distributions of individual parameters of the 2002 H II regions constituting one of the composite star-forming galaxies built as described in Section 4.1 (the meaning of the different parameters is described in Table 3). The corresponding ‘baseline’ parameters are shown as vertical black dashed lines. The vertical red lines show the average, Hα-weighted H II-region parameters computed as described in Section 5.1. The dashed yellow line in the $A_V$ panel shows the value of this parameter recovered from the global Hα/Hβ ratio, as described in Section 4.2.</figcaption></figure>

and the fact that the rate of H-ionizing photons, $Q_0$, scales as $U^3$ (Section 2.1), implies a deficit of galaxies with low O/H. Similarly, the geometric form factor $f_r$ also depends on $\log U$ through $Q_0$ and the limit on LHα. The shapes of the X/O distributions in Fig. 8 originate from the dependence of $(X/H)^{\text{tot}}$ on $(O/H)^{\text{tot}}$ according to Nicholls et al. (2017). Most notably, the sharp increase in $(C/O)^{\text{tot}}$ with $(O/H)^{\text{tot}}$ over the range $8.3 \lesssim 12 + \log(O/H) \lesssim 9.0$ and the plateau values outside this range imply that for $\log(O/H) \lesssim -3.7$ (corresponding to most of the range in Fig. 8), galaxies stack up at low C/O and are more sparsely sampled at higher C/O, corresponding to $-3.7 \lesssim 12 + \log(O/H) \lesssim 3.0$. The trend is less

marked for N/O, due to the adoption of a smooth dependence of $(N/O)^{\text{tot}}$ on $(O/H)^{\text{tot}}$ by Nicholls et al. (2017, their equation 9), while $(S/O)^{\text{tot}}$, $(Ne/O)^{\text{tot}}$, and $(Ar/O)^{\text{tot}}$ are assumed to vary at most negligibly with $12 + \log(O/H)$ (their table 2). Other parameters, such as $n_H$ and $A_V$, have more symmetric distributions, sharply peaked near the middle of the fixed ranges over which individual H II-region parameters are drawn for all galaxies. We note the uneven distribution of $N_{HII}$, which arises from the random drawing of LHα combined with the imposed Hα luminosity function (equation 6).

It is of particular interest to check how the 250 000 composite star-forming galaxies built in this way populate standard line-
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Distributions of various properties of the 250 000 composite star-forming galaxies built as described in Section 4 and Table 3. Each galaxy contains between $N_{HII} = 100$ and about 3000 constituent H II regions, from which average H II-region parameters are computed by weighting the individual H II-region parameters by the intrinsic Hα luminosities of these regions (see Section 5.1 for more details and comments on the distribution shapes).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Distribution of the 250 000 composite star-forming galaxies of Fig. 8 in various line-diagnostic diagrams, colour-coded according to oxygen abundance as indicated on the right-hand side. In the [O III] λ5007/Hβ-versus-[N II] λ6584/Hα (BPT) diagram, the separation criteria between AGNs and star-forming galaxies by Kewley et al. (2001, solid line), Kauffmann et al. (2003, dotted line), and Stasińska et al. (2006, dashed line) are shown. Similar criteria by Kewley et al. (2001) appear as solid lines in the [O III] λ5007/Hβ-versus-[S II] λ6716, 31/Hα and [O III] λ5007/Hβ-versus-[O I] λ6300/Hα (VO87) diagrams, and by Backhaus et al. (2022) in the [O III] λ5007/Hβ-versus-[Ne III] λ3869/[O II] λ3727 (Zeimann et al. 2015) diagram (see Section 5.1 for details).</figcaption></figure>

diagnostic diagrams. This is shown in Fig. 9 for different diagrams involving [O II] λ3727, [Ne III] λ3869, Hβ, [O III] λ5007, [O I] λ6300, Hα, [N II] λ6584, [S II] λ6725, P₁, and P₂ (footnote 5), inspired from various previous studies (BPT; VO87; Zeimann et al. 2015; Ji & Yan 2022). In the [O III] λ5007/Hβ-versus-[N II] λ6584/Hα (BPT) diagram, we show the criteria of Kewley et al. (2001, solid line), Kauffmann et al. (2003, dotted line), and Stasińska et al. (2006, dashed line) to separate active galactic nucleus (AGN)-dominated from star-forming galaxies. Solid lines show similar criteria by Kewley et al. (2001) in the [O III] λ5007/Hβ-versus-[S II] λ6725/Hα and [O III] λ5007/Hβ-versus-[O I] λ6300/Hα (VO87) diagrams, and by Backhaus et al. (2022) in the [O III] λ5007/Hβ-versus-[Ne III] λ3869/[O II] λ3727 (Zeimann et al. 2015) diagram.

The composite star-forming galaxies' models fall roughly on the expected side of these criteria in all diagrams, except in the case of the [O III] λ5007/Hβ-versus-[Ne III] λ3869/[O II] λ3727 diagram. This is not surprising, as the Backhaus et al. (2022) criterion in this diagram marks the region below which AGN-dominated galaxies are unlikely to be found, but does not exclude star-forming galaxies on either side (see their fig. 6 and section 3.2, as well as Trump et al. 2023; Hu et al. 2024).

The models in Fig. 9, which are colour-coded according to oxygen abundance, show well-defined sequences in 12 + log (O/H) in all diagrams, as expected from many previous studies (e.g. Pagel et al. 1979; McGaugh 1991; Charlot & Longhetti 2001; Kewley et al. 2001). Fig. D1 shows analogues to Fig. 9 colour-coded according to other quantities: the DIG contribution to Hβ luminosity (DIGFrac), and the offsets $\rho(\delta_{N/O})$ and $\rho(\delta_U)$ of the baseline galactic N/O and log U with respect to canonical relations (Section 4.1). These show, for example, that star-forming galaxies located on the AGN side of the Kewley et al. (2001) relation in the BPT diagram tend to have large DIGFrac and $\rho(\delta_{N/O})$.

## 5.2 O/H determination from integrated line ratios

We now examine how the oxygen abundance of a composite star-forming galaxy derived using the classical 'direct method' (e.g. Aller & Liller 1959; Peimbert 1967; Peimbert & Costero 1969) differs from the true Hα-weighted oxygen abundance of its constituent H II regions.

<table><caption>**Table 4.** Atomic data sets used for collisionally excited lines in the CLOUDY and PYNEB calculations presented in this paper.</caption><thead><tr><th>Ion</th><th>Transition probabilities</th><th>Collision strengths</th></tr></thead><tbody><tr><td>N<sup>+</sup></td><td>Tachiev & Froese Fischer (2001)</td><td>Tayal (2011)</td></tr><tr><td>O<sup>+</sup></td><td>Froese Fischer & Tachiev (2004)</td><td>Kisielius et al. (2009)</td></tr><tr><td>O<sup>2+</sup></td><td>Tachiev & Froese Fischer (2001)</td><td>Storey, Sochi & Badnell (2014)</td></tr><tr><td>S<sup>+</sup></td><td>Kisielius et al. (2014)</td><td>Tayal & Zatsarinny (2010)</td></tr></tbody></table>

We derive the electronic temperature and density of the emitting gas in a given galaxy using the PYNEB code (version 1.1.19; Luridiana, Morisset & Shaw 2015; Morisset et al. 2020), adopting the same atomic data as used to generate the CLOUDY v17.03 photoionization models in Section 2.2, whose references are listed in Table 4. PYNEB incorporates the effects of collisional excitation for lines emitted by metals. We use here a simplified two-zone ionization structure. The electronic temperatures of the N⁺ and O²⁺ zones are derived from the ratios of auroral to nebular forbidden lines [N II] λ5755/λ6584 and [O III] λ4363/λ5007, and the electronic density of the S⁺ zone from the ratio of forbidden lines [S II] λ6716/λ6731. PYNEB computes the O⁺/H⁺ and O²⁺/H⁺ ionic abundances assuming that the electronic density is the same everywhere and that the electronic temperature of the N⁺ zone is representative of that of the O⁺ zone. The sum of O⁺/H⁺ and O²⁺/H⁺ gives the total O/H abundance of the galaxy, assuming that the contributions of O⁰ and O³⁺ are negligible. No correction is applied to take into account emission process other than collisional excitation (e.g. recombination or fluorescence).

In the upper panel of Fig. 10, we show a density map of the logarithmic difference between O/H derived in this way and the true, Hα-weighted O/H as a function of 12 + log (O/H), for the subset of 207 750 composite star-forming galaxies of Fig. 8 with significant [O III] λ4363 emission, i.e. with log([O III] λ4363/Hβ) > −3. We notice that at the lowest metallicities, i.e. for 12 + log (O/H) ≲ 7.5, the direct method tends to slightly underestimate the oxygen abundance, by about 0.05 dex. This result, reminiscent of the early finding
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Distribution of the logarithmic difference between O/H derived from the direct method and true, Hα-weighted O/H as a function of 12 + log (O/H), for the subset of 207 750 composite star-forming galaxies of Fig. 8 with log([O III] λ4363/Hβ) > −3. Top: Colour-coded as a function of model number density. Bottom: Colour-coded as a function of root-mean-square t of H II-region temperature fluctuations (equation 7).</figcaption></figure>

by Kobulnicky et al. (1999), is likely caused by the neglect of O³⁺ in our PYNEB estimates, which starts to be significant as the ionization parameter increases with decreasing O/H (equation 3).

Most notably, Fig. 10 shows that direct-method determinations of O/H in metal-rich galaxies can seriously underestimate the true O/H, by factors of up to several. Several origins for this phenomenon have been proposed, among them the neglect by the direct method, at high metallicity (i.e. low electronic temperature), of the contribution of recombination to the [N II] λ5755 and [O III] λ4363 auroral lines (Gómez-Llanos et al. 2020, and references therein) and the fact that any temperature gradient within the nebula preferentially enhances [O III] λ5007 emission in low-temperature regions, while [O III] λ4363 is more prominent in high-temperature regions (Stasińska 2005). Moreover, our modelling of the electronic temperature (T_e) distribution of H II regions in composite galaxies offers a unique opportunity to study the influence of T_e fluctuations on O/H estimates of star-forming galaxies using the direct method (even though our CLOUDY calculations do not include such fluctuations within individual regions). By analogy with Peimbert (1967)'s parametrization of temperature fluctuations inside H II regions, we define the root-mean-square t of temperature fluctuations (weighted by Hα luminosity) among H II regions in a composite star-forming galaxy by

$$t^2 = \left\langle \left( \frac{T_e - T_0}{T_0} \right)^2 \right\rangle_{\text{LH}\alpha}, \qquad (7)$$

where LHα is the Hα luminosity of the region of electronic temperature T_e, and

$$T_0 = \langle T_e \rangle_{\text{LH}\alpha} \qquad (8)$$

is the mean, Hα-weighted temperature of H II regions in the galaxy (Section 5.1). The t² distribution of the galaxies in our sample is shown in the bottom-right panel of Fig. 8.

The lower panel of Fig. 10 shows the analogue of the upper panel, but colour-coded according to log(t²). Remarkably, this plot reveals that the galaxies for which the direct method most strongly underestimates O/H at high metallicities exhibit the strongest T_e fluctuations, with deficits of over 0.2 dex corresponding to t² > 0.03. This result suggests that the spread in T_e among the constituent H II regions of a galaxy might be an important factor causing the direct method to underestimate O/H. This echoes the conclusions drawn by Cameron et al. (2023) from the radiation-hydrodynamics simulation of an isolated dwarf galaxy.

Hence, Fig. 10 illustrates the uncertainties and biases that can affect O/H determinations when obtained for observations in which light from numerous individual H II regions and DIG are collected together. A detailed investigation of how these effects influence abundance determinations in samples of observed star-forming galaxies (and, a fortiori, abundances derived from stacking several galaxy spectra), as well as the impact of noise on these results, will be addressed in a forthcoming study.

## 6 SUMMARY AND PERSPECTIVES

We have developed a novel approach to model the nebular emission from star-forming galaxies, incorporating the contributions of numerous individual H II regions and a DIG component (considering only DIG produced by HOLMES and excluding contributions from photon leakage from H II regions, which may dominate in at least some star-forming galaxies; see Lugo-Aranda et al. 2024). The approach leverages an ML-based regression ANN, trained on a grid of models generated by the photoionization code CLOUDY, to predict the emission-line properties of H II regions under various physical conditions.

The ANN can perform 1 million predictions of 14 line ratios, electron temperature and Hα luminosity in less than 0.2 s from a sampling of the 13 input parameters, on an NVIDIA RTX A5000 GPU graphics card. This CLOUDY substitute can be applied in scenarios requiring the generation of numerous models, such as MCMC simulations and generic algorithms involving multiple model iterations.

We constructed 250 000 synthetic star-forming galaxies by summing, for each galaxy, the emission from between 100 and about 3000 H II regions following an Hα luminosity function of slope α_LF = −1.9 (Rousseau-Nepton et al. 2018; Santoro et al. 2022), each with different physical parameters and dust attenuation, and a component of diffuse gas photoionized by HOLMES (Flores-Fajardo et al. 2011). The distributions of H II-region parameters are tunable and chosen to reproduce observed trends, the relationships between O/H, log U, C/O, N/O, Ne/O, S/O, and Ar/O being considered as adjustable parameters changing from one galaxy to another. The integrated properties for each synthetic galaxy are defined as Hα-weighted mean values of the corresponding properties of the contributing H II regions. The integrated properties of the resulting 250 000 galaxies lie roughly in the regions expected to be populated by star-forming galaxies in standard line-ratio diagnostic diagrams involving [O II] λ3727, [Ne III] λ3869, Hβ, [O III] λ5007, [O I] λ6300, Hα, [N II] λ6584, and [S II] λ6725 (BPT; VO87; Zeimann et al. 2015;
Ji & Yan 2022). This approach enables one to explore how different parameter distributions affect the integrated emission-line properties of composite star-forming galaxies.

As a practical application, we examined how the oxygen abundances of composite star-forming galaxies determined using the traditional 'direct method' (e.g. Aller & Liller 1959; Peimbert 1967; Peimbert & Costero 1969) compare to the actual Hα-weighted abundances derived from the properties of the individual H II regions that constitute galaxies. We find that the direct method exhibits systematic biases at both low [12 + log (O/H) < 7.5] and high [12 + log (O/H) > 8.5] metallicities, with a more pronounced effect at the higher end, resulting in underestimates of O/H in both cases. The bias at low metallicity (reminiscent of that found in the pioneer study of Kobulnicky et al. 1999) arises from the common omission of O³⁺ in O/H estimates achieved in the two-zone (O⁺/O²⁺) approximation using the standard PyNeb package (Luridiana et al. 2015; Morisset et al. 2020). At high metallicities, the bias is most severe in galaxies with the largest spread in electronic temperatures among their constituent H II regions, aligning with the findings of Cameron et al. (2023) from radiation-hydrodynamics simulations of an isolated dwarf galaxy. These findings highlight the need to account for the composite nature of star-forming galaxies when deriving physical parameters from the observed integrated nebular emission.

The results presented in this paper are primarily illustrative, the main objective being to introduce our novel approach to modelling efficiently the integrated properties of composite star-forming galaxies. We have limited here our investigation to the optical emission-line properties of galaxies, but similar analyses could be extended to include ultraviolet and infrared emission-line diagnostics (e.g. Cormier et al. 2019; Hirschmann et al. 2019). This would involve training a new ANN to predict the relevant line luminosities and generating composite galaxy models by summing these additional line emissions among the constituent H II regions.

In a forthcoming paper, we will apply the approach introduced here to fitting observed line ratios of star-forming galaxies. We plan to explore a new parametrization of the ANN model to avoid generating unrealistic combinations of input parameters, which lead to improbable Hα emissions. We also plan to use more efficient distributions for the input parameters, such as the Latin Hypercube, the Sobol, or the Halton ones. Additionally, we may update the training using models obtained from the CLOUDY C23.01 version.

## ACKNOWLEDGEMENTS

The authors would like to thank the anonymous referee for their valuable report, which significantly contributed to improving the quality of this paper. The initial stages of this work were conceived and developed during a CM sabbatical at the Institut d'Astrophysique de Paris (IAP) to collaborate with SC. CM gratefully acknowledges the PASPA grant that supported this sabbatical. The authors thank Grazyna Stasińska for valuable discussions. This work used v2.2.1 of the BPASS models as described in Eldridge et al. (2017) and Stanway & Eldridge (2018). CM acknowledges the support of UNAM/DGAPA/PAPIIT grants IN101220 and IG101223.

**Software:** In this study, the following PYTHON libraries are used: NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), MATPLOTLIB (Hunter 2007), ASTROPY (The Astropy Collaboration 2022), PYNEB (Luridiana et al. 2015; Morisset et al. 2020), PYCLOUDY (Morisset et al. 2014), SQLALCHEMY (Bayer 2012), PANDAS (McKinney 2010; The pandas development team 2020), SCIKIT-LEARN (Pedregosa et al. 2011), and TENSORFLOW (Abadi et al. 2015). This paper has been edited using Overleaf facilities (https://www.overleaf.com/).

## DATA AVAILABILITY

The PYTHON scripts and notebooks used for this research are available from the GitHub repository at https://github.com/Morisset/pysyngal (locked until publication of this and the forthcoming paper II). The repository offers a downloadable version of the trained ANN, as well as the training set and training script. The original CLOUDY models can be found on 3MDB_17 under reference HII_24.

## REFERENCES

* Abadi M. et al., 2015, preprint (arXiv:1605.08695)

* Aller L. H., Liller W., 1959, ApJ, 130, 45

* Alloin D., Collin-Souffrin S., Joly M., Vigroux L., 1979, A&A, 78, 200

* Astropy Collaboration, 2022, ApJ, 935, 167

* Backhaus B. E. et al., 2022, ApJ, 926, 161

* Badenas-Agusti M., Viaña J., Vanderburg A., Blouin S., Dufour P., Xu S., Sha L., 2024, MNRAS, 529, 1688

* Baldwin A., Phillips M. M., Terlevich R., 1981, PASP, 93, 817

* Baron D., 2019, preprint (arXiv:1904.07248)

* Bayer M., 2012, in Brown A., Wilson G., eds, The Architecture of Open Source Applications, Vol. II: Structure, Scale, and a Few More Fearless Hacks. Available at: aosabook.org

* Bird S., Rogers K. K., Peiris H. V., Verde L., Font-Ribera A., Pontzen A., 2019, J. Cosmol. Astropart. Phys., 2019, 050

* Bron E. et al., 2021, A&A, 645, A28

* Byler N., Dalcanton J. J., Conroy C., Johnson B. D., Levesque E. M., Berg D. A., 2018, ApJ, 863, 14

* Byrne C. M., Stanway E. R., Eldridge J. J., McSwiney L., Townsend O. T., 2022, MNRAS, 512, 5329

* Cameron A. J., Katz H., Rey M. P., 2023, MNRAS, 522, L89

* Carton D. et al., 2017, MNRAS, 468, 2140

* Chabrier G., 2003, ApJ, 586, L133

* Charlot S., Longhetti M., 2001, MNRAS, 323, 887

* Cormier D. et al., 2019, A&A, 626, A23

* DeRose J., Chen S.-F., White M., Kokron N., 2022, J. Cosmol. Astropart. Phys., 2022, 056

* Dopita M. A., Sutherland R. S., Nicholls D. C., Kewley L. J., Vogt F. P. A., 2013, ApJS, 208, 10

* Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Aust., 34, e058

* Ferguson J. W., Korista K. T., Baldwin J. A., Ferland G. J., 1997, ApJ, 487, 122

* Ferland G. J. et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385

* Fernández V., Amorín R., Sanchez-Janssen R., del Valle-Espinosa M. G., Papaderos P., 2023, MNRAS, 520, 3576

* Fitzpatrick E. L., 1999, PASP, 111, 63

* Flores-Fajardo N., Morisset C., Stasińska G., Binette L., 2011, MNRAS, 415, 2182

* Froese Fischer C., Tachiev G., 2004, At. Data Nucl. Data Tables, 87, 1

* Galliano F., 2018, MNRAS, 476, 1445

* Gómez-Llanos V., Morisset C., García-Rojas J., Jones D., Wesson R., Corradi R. L. M., Boffin H. M. J., 2020, MNRAS, 498, L82

* Grassi T., Nauman F., Ramsey J. P., Bovino S., Picogna G., Ercolano B., 2022, A&A, 668, A139

* Gutkin J., Charlot S., Bruzual G., 2016, MNRAS, 462, 1757

* Harris C. R. et al., 2020, Nature, 585, 357

* Heckman T. M. et al., 2011, ApJ, 730, 5

* Hendrycks D., Gimpel K., 2016, preprint (arXiv:1606.08415)

* Hirschmann M., Charlot S., Feltre A., Naab T., Somerville R. S., Choi E., 2019, MNRAS, 487, 333

* Ho I. T., 2019, MNRAS, 485, 3569

* Holdship J., Viti S., Haworth T. J., Ilee J. D., 2021, A&A, 653, A76

* Hu W. et al., 2024, ApJ, 971, 21

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

* Jensen E. B., Strom K. M., Strom S. E., 1976, ApJ, 209, 748

* Ji X., Yan R., 2022, A&A, 659, A112
Katz H. et al., 2023, MNRAS, 518, 592

Kauffmann G. et al., 2003, MNRAS, 346, 1055

Kewley L. J., Dopita M. A., 2002, ApJS, 142, 35

Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183

Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J., 2001, ApJ, 556, 121

Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511

Kingma D. P., Ba J., 2014, preprint (arXiv:1412.6980)

Kisielius R., Storey P. J., Ferland G. J., Keenan F. P., 2009, MNRAS, 397, 903

Kisielius R., Kulkarni V. P., Ferland G. J., Bogdanovich P., Lykins M. L., 2014, ApJ, 780, 76

Kobulnicky H. A., Kennicutt R. C., Jr, Pizagno J. L., 1999, ApJ, 514, 544

Korobov A., 1959, Dokl. Akad. Nauk SSSR, 124, 1207

Lebouteiller V. et al., 2024, preprint (arXiv:2412.15860)

Lecroq M. et al., 2024, MNRAS, 527, 9480

Li Y., Leja J., Johnson B. D., Tacchella S., Davies R., Belli S., Park M., Emami R., 2024, preprint (arXiv:2405.04598)

Liu D., Nocedal J., 1989, Math. Program., 45, 503

Lugo-Aranda A. Z., Sánchez S. F., Barrera-Ballesteros J. K., López-Cobá C., Espinosa-Ponce C., Galbany L., Anderson J. P., 2024, MNRAS, 528, 6099

Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42

McClymont W. et al., 2024, MNRAS, 532, 2016

McGaugh S. S., 1991, ApJ, 380, 140

McKinney W., 2010, in van der Walt S., Millman J., eds, Proc. 9th Python in Science Conference, p. 56

Maltsev K., Schneider F. R. N., Röpke F. K., Jordan A. I., Qadir G. A., Kerzendorf W. E., Riedmiller K., van der Smagt P., 2024, A&A, 681, A86

Marino R. A. et al., 2013, A&A, 559, A114

Martínez-Paredes M., Bruzual G., Morisset C., Kim M., Meléndez M., Binette L., 2023, MNRAS, 525, 2916

Mast D. et al., 2014, A&A, 561, A129

Mootoovaloo A., Jaffe A. H., Heavens A. F., Leclercq F., 2022, Astron. Comput., 38, 100508

Morisset C., 2014, in Morisset C., Delgado-Inglada G., Torres-Peimbert S., eds, Proc. Asymmetrical Planetary Nebulae VI Conference. Mexico

Morisset C., Delgado-Inglada G., Flores-Fajardo N., 2015, Rev. Mex. Astron. Astrofis., 51, 103

Morisset C. et al., 2016, A&A, 594, A37

Morisset C., Luridiana V., García-Rojas J., Gómez-Llanos V., Bautista M., Mendoza C., 2020, Atoms, 8, 66

Nicholls D. C., Sutherland R. S., Dopita M. A., Kewley L. J., Groves B. A., 2017, MNRAS, 466, 4403

Osterbrock D. E., Ferland G. J., 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei. University Science Books, Sausalito, CA

Pagel B. E. J., Edmunds M. G., Blackwell D. E., Chun M. S., Smith G., 1979, MNRAS, 189, 95

Palud P. et al., 2023, A&A, 678, A198

Pedregosa F. et al., 2011, J. Mach. Learn. Res., 12, 2825

Peimbert M., 1967, ApJ, 150, 825

Peimbert M., Costero R., 1969, Bol. Obs. Tonantzintla Tacubaya, 5, 3

Pérez-Montero E., 2017, PASP, 129, 043001

Pérez-Montero E., Amorín R., Sánchez Almeida J., Vílchez J. M., García-Benito R., Kehrig C., 2021, MNRAS, 504, 1237

Postnikova V. K., Bizyaev D., 2023, Astron. Lett., 49, 151

Ramambason L. et al., 2022, A&A, 667, A35

Ramirez S., Icaza-Lizaola M., Fromenteau S., Vargas-Magaña M., Aviles A., 2024, J. Cosmol. Astropart. Phys., 2024, 049

Rautio R. P. V., Watkins A. E., Comerón S., Salo H., Díaz-García S., Janz J., 2022, A&A, 659, A153

Reza M., Zhang Y., Nord B., Poh J., Ciprijanovic A., Strigari L., 2022, in ICML 2022 Workshop on Machine Learning for Astrophysics. p. 20

Richardson C. T., Allen J. T., Baldwin J. A., Hewett P. C., Ferland G. J., Crider A., Meskhidze H., 2016, MNRAS, 458, 988

Rousseau-Nepton L., Robert C., Martin R. P., Drissen L., Martin T., 2018, MNRAS, 477, 4152

Sanders R. L., Shapley A. E., Zhang K., Yan R., 2017, ApJ, 850, 136

Santoro F. et al., 2022, A&A, 658, A188

Schaefer A. L., Tremonti C., Belfiore F., Pace Z., Bershady M. A., Andrews B. H., Drory N., 2020, ApJ, 890, L3

Schaefer D., Fragos T., Izotov Y. I., 2019, A&A, 622, L10

Senchyna P., Stark D. P., Mirocha J., Reines A. E., Charlot S., Jones T., Mulchaey J. S., 2020, MNRAS, 494, 941

Smirnov-Pinchukov G. V., Molyarova T., Semenov D. A., Akimkin V. V., van Terwisga S., Francheschi R., Henning T., 2022, A&A, 666, L8

Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75

Stasińska G., 2005, A&A, 434, 507

Stasinska G., 2019, preprint (arXiv:1906.04520)

Stasińska G., Cid Fernandes R., Mateus A., Sodré L., Asari N. V., 2006, MNRAS, 371, 972

Stasińska G., Izotov Y., Morisset C., Guseva N., 2015, A&A, 576, A83

Storey P. J., Sochi T., Badnell N. R., 2014, MNRAS, 441, 3028

Tachiev G., Froese Fischer C., 2001, Can. J. Phys., 79, 955

Tayal S. S., 2011, ApJS, 195, 12

Tayal S. S., Zatsarinny O., 2010, ApJS, 188, 32

The pandas, 2020, pandas-dev/pandas: Pandas, Zenodo, available at: https://zenodo.org/records/13819579

Trump J. R. et al., 2023, ApJ, 945, 35

Umeda H., Ouchi M., Nakajima K., Isobe Y., Aoyama S., Harikane Y., Ono Y., Matsumoto A., 2022, ApJ, 930, 37

Vale Asari N., Stasińska G., Morisset C., Cid Fernandes R., 2016, MNRAS, 460, 1739

Veilleux S., Osterbrock D. E., 1987, ApJS, 63, 295

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Wu J. F., Boada S., 2019, MNRAS, 484, 4683

Zeimann G. R. et al., 2015, ApJ, 798, 29

## APPENDIX A: USING A NEURAL NETWORK AS A REGRESSOR FOR A GRID OF CLOUDY H II REGION MODELS

The ML method used to substitute CLOUDY in the prediction of emission-line ratios is of ANN type. This ANN is implemented using the PYTHONTENSORFLOW library (http://tensorflow.org/). The neural architecture counts three hidden layers of 256, 512, and 256 activation cells, respectively. All activation functions are of GELU type (Hendrycks & Gimpel 2016), which combines the fast-convergence benefits of ReLU (Rectified Linear Unit) and ELU (Exponential Linear Unit) functions with the regularization effect of a dropout, achieved by randomly setting some activations to zero during training. The GELU function is infinitely differentiable. The input layer accepts a 13-dimensional tensor, and the output layer generates a 16-dimensional tensor (see Section 3 for the physical parameters and observables associated with the input and output). The loss function is the negative logarithmic likelihood of the Gaussian distribution where the mean and standard deviation are estimated. This choice allows us to describe the data as noisy observations with heteroscedastic Gaussian noise. In the homoscedastic case, only the mean is estimated, and one recovers the classical Mean Squared Error (MSE) loss.

The optimizer is Adam (Kingma & Ba 2014) with parameters $\beta_1 = 0.9$ and $\beta_2 = 0.999$. A learning-rate scheduler is tailored using an alternation of exponential decays, starting at epoch 500. The goal is to maximize the rate at which the loss decreases while ensuring loss minimization with a learning rate ending at almost zero. Its value starts at 0.001 and ends at approximately $10^{-7}$.

Training is performed with a batch size of 2000, using 80 per cent of the grid of 1 million CLOUDY models described in Section 2, with ionizing-cluster age $t_c < 7$ Myr and $H\beta$ Frac > 0.65. We explore the impact of changing the training-set size, batch size, and number of epochs. The results are presented in Table A1 for seven combinations of training-set sizes and number of epochs. The standard deviations
<table><caption>**Table A1.** Performance of different ANNs with different training-set sizes and numbers of epochs. The training set consists of 80 per cent of the total data, while 20 per cent is used to evaluate the fit quality. The table lists the standard deviation (times 100) of the difference between predicted and true values (both being logarithm values of line ratios, $T_e$ or $H\alpha$) for each element of the prediction vector in the test set. The last row shows the mean of these 16 standard deviations.</caption><thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>Total set size (×10<sup>5</sup>)</td><td>2</td><td>4</td><td>6</td><td>6</td><td>8</td><td>10</td><td>10</td><td>10</td></tr><tr><td>Epochs (×10<sup>3</sup>)</td><td>2</td><td>2</td><td>2</td><td>10</td><td>2</td><td>2</td><td>8</td><td>8</td></tr><tr><td>Batch size (×10<sup>3</sup>)</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>0.5</td></tr><tr><td>Training time (mn)</td><td>5.8</td><td>11.6</td><td>14.5</td><td>72.1</td><td>19.4</td><td>23.5</td><td>88.3</td><td>345.0</td></tr><tr><td>ANN output</td><td></td><td></td><td></td><td colspan="4">Standard deviation (×100)</td><td></td></tr><tr><td>Hα/Hβ</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.1</td><td>0.2</td><td>0.2</td><td>0.1</td><td>0.1</td></tr><tr><td>[O III] λ5007/Hβ</td><td>4.7</td><td>3.4</td><td>2.6</td><td>2.2</td><td>2.1</td><td>1.8</td><td>1.3</td><td>1.3</td></tr><tr><td>[N II] λ6584/Hβ</td><td>4.4</td><td>3.2</td><td>2.5</td><td>2.2</td><td>2.2</td><td>1.9</td><td>1.6</td><td>1.5</td></tr><tr><td>[O II] λ3727/Hβ</td><td>5.0</td><td>3.8</td><td>3.1</td><td>2.8</td><td>2.6</td><td>2.2</td><td>1.7</td><td>1.7</td></tr><tr><td>[O III] λ4363/λ5007</td><td>3.6</td><td>2.9</td><td>2.3</td><td>2.1</td><td>2.0</td><td>1.8</td><td>1.3</td><td>1.3</td></tr><tr><td>[N II] λ5755/λ6584</td><td>3.4</td><td>2.6</td><td>2.1</td><td>1.9</td><td>1.8</td><td>1.6</td><td>1.4</td><td>1.3</td></tr><tr><td>[S II] λ6725/Hβ</td><td>4.4</td><td>3.1</td><td>2.5</td><td>2.3</td><td>2.2</td><td>1.9</td><td>1.7</td><td>1.5</td></tr><tr><td>[S III] λ9530/Hβ</td><td>3.1</td><td>2.2</td><td>1.7</td><td>1.5</td><td>1.4</td><td>1.2</td><td>0.9</td><td>0.9</td></tr><tr><td>[Ar III] λ7135/Hβ</td><td>3.6</td><td>2.6</td><td>2.0</td><td>1.7</td><td>1.7</td><td>1.4</td><td>1.1</td><td>1.0</td></tr><tr><td>[Ne III] λ3868/Hβ</td><td>6.2</td><td>4.4</td><td>3.4</td><td>2.9</td><td>2.8</td><td>2.4</td><td>1.7</td><td>1.7</td></tr><tr><td>[S II] λ6716/λ6731</td><td>0.5</td><td>0.4</td><td>0.4</td><td>0.2</td><td>0.3</td><td>0.3</td><td>0.2</td><td>0.2</td></tr><tr><td>He I λ5876/Hβ</td><td>1.1</td><td>1.0</td><td>0.9</td><td>0.8</td><td>0.9</td><td>0.8</td><td>0.8</td><td>0.8</td></tr><tr><td>[S III] λ6312/λ9530</td><td>3.4</td><td>2.5</td><td>2.0</td><td>1.7</td><td>1.6</td><td>1.4</td><td>1.1</td><td>1.1</td></tr><tr><td>[O II] λ7325/Hβ</td><td>5.5</td><td>4.2</td><td>3.4</td><td>3.1</td><td>2.8</td><td>2.3</td><td>1.8</td><td>1.8</td></tr><tr><td>T<sub>e</sub></td><td>1.0</td><td>0.8</td><td>0.6</td><td>0.5</td><td>0.6</td><td>0.5</td><td>0.4</td><td>0.4</td></tr><tr><td>LHα</td><td>2.1</td><td>1.9</td><td>1.4</td><td>0.9</td><td>1.4</td><td>1.4</td><td>0.8</td><td>0.8</td></tr><tr><td>Mean</td><td>3.3</td><td>2.5</td><td>1.9</td><td>1.7</td><td>1.7</td><td>1.5</td><td>1.1</td><td>1.1</td></tr></tbody></table>

listed in the table are obtained for test points in the triangle defined in Section 3.2.

Training times are provided for an NVIDIA RTX A5000 graphics card. Standard deviations of the difference between predicted and true values for each output tensor element are reported. The quality of predictions improves with larger training sets and more epochs. The ANN trained with 800 000 models over 8000 epochs achieves an impressive prediction capacity of close to 1 per cent in dex. This ANN is used for the analysis in this paper. The input tensor is scaled using the SCIKIT-LEARN StandardScaler function. The ANN is managed through the AI4NEB library.

## APPENDIX B: IMPACT OF SELECTED PARAMETERS ON EMISSION-LINE RATIOS OF H II-REGION MODELS

Figs B1 and B2 illustrate how variations in ionizing-cluster age, gas density, and the HβFrac and $\xi_d$ parameters (from top to bottom) impact the emission-line ratios in the same diagrams as in the six rightmost panels of Fig. 5. These figures reveal the complexity of these effects and the potential degeneracies that arise when attempting to reproduce the observed line ratios.
<figure><img src="image_13.png" /><figcaption>**Figure B1.** Same as the six rightmost panels of Fig. 5, but colour-coded according to ionizing-cluster age $t_c$ (top) and gas density $n_H$ (bottom).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure B2.** Same as the six rightmost panels of Fig. 5, but colour-coded according to Hβ-luminosity ratio HβFrac of truncated versus radiation-bounded model (top) and dust-to-metal mass ratio ξd (bottom).</figcaption></figure>
## APPENDIX C: USING A NEURAL NETWORK AS A REGRESSOR FOR A GRID OF CLOUDY DIG MODELS

The DIG is assumed to primarily consist of regions photoionized by HOLMES. To represent the DIG, we use the photoionization models of Martínez-Paredes et al. (2023), obtained from the 3MdB_17 data base under the reference 'CB_19'. From this data set, we select 9000 models with ages exceeding 1 Myr and H densities $n_H \le 10 \text{ cm}^{-3}$. We randomly split these models into a training set (80 per cent) and a test set (20 per cent).

We train an ANN, referred to as 'DIG-ANN', using six input parameters: the age of the ionizing stars, drawn uniformly between 1 and 10 Gyr; the O/H, C/O, and N/O abundances in the DIG, assumed to match the global abundances of the composite galaxy (Section 4.1); the ionization parameter, drawn uniformly in the range $-4 \le \log U \le -3$; and the quantity HβFrac, drawn uniformly between 0.4 and 1.0. The outputs are the 14 emission-line ratios described in Section 3.

The DIG-ANN is implemented using the SCIKIT-LEARN library (Pedregosa et al. 2011), with three hidden layers comprising 50, 100, and 50 tanh-based perceptrons, respectively. Convergence is

achieved using the LBFGS (Limited memory Broyden-Fletcher-Goldfarb-Shanno) solver (Liu & Nocedal 1989). The DIG-ANN is managed through the AI4NEB library (footnote 4).

The standard deviations of the differences between the predicted and the test logarithmic values of the emission-line ratios are of the order of 0.02, comparable to the performance of the main ANN reported in Table A1.

## APPENDIX D: COMPOSITE STAR-FORMING GALAXIES IN EMISSION-LINE DIAGRAMS COLOUR-CODED BY DIFFERENT PHYSICAL PARAMETERS

Fig. D1 shows analogues of Fig. 9 displaying the distribution of the 250 000 composite star-forming galaxies in various line-ratio diagnostic diagrams, but colour-coded according to the DIG contribution to Hβ luminosity (DIGFrac) and the offsets $\rho(\delta_{N/O})$ and $\rho(\delta_U)$ of the baseline galactic N/O and $\log U$ relative to canonical relations (Section 4.1). The results indicate, for example, that star-forming galaxies on the AGN side of the Kewley et al. (2001) relation in the [O III] $\lambda5007/H\beta$-versus-[N II] $\lambda6584/H\alpha$ (BPT) diagram generally exhibit high DIGFrac and $\rho(\delta_{N/O})$.
<figure><img src="image_15.png" /><figcaption>**Figure D1.** Same as Fig. 9, but colour-coded according to DIG contribution to Hβ luminosity (DIGFrac; top panels) and offsets $\rho(\delta_{N/O})$ (middle panels) and $\rho(\delta_{U})$ (bottom panels) of the baseline galactic N/O and log U relative to the canonical relations described in Section 4.1.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.