{'fraction_non_alphanumeric': 0.02743505249444402, 'fraction_numerical': 0.008791040364341, 'mean_word_length': 4.615639985245297, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
WARPFIELD-EMP: The Self-Consistent Prediction of Emission Lines from Evolving HII Regions in Dense Molecular Clouds

We present the {\sc warpfield} emission predictor, {\sc warpfield-emp}, which couples the 1D stellar feedback code {\sc warpfield} with the {\sc cloudy} \hii region/PDR code and the {\sc polaris} radiative transfer code, in order to make detailed predictions for the time-dependent line and continuum emission arising from the H{\sc ii} region and PDR surrounding an evolving star cluster. {\sc warpfield-emp} accounts for a wide range of physical processes (stellar winds, supernovae, radiation pressure, gravity, thermal conduction, radiative cooling, dust extinction etc.) and yet runs quickly enough to allow us to explore broad ranges of different cloud parameters. We compare the results of an extensive set of models with SITELLE observations of a large sample of \hii regions in NGC~628 and find very good agreement, particularly for the highest signal-to-noise observations. We show that our approach of modeling individual clouds from first principles (instead of in terms of dimensionless quantities such as the ionization parameter) allows us to avoid long-standing degeneracies in the interpretation of \hii region diagnostics and enables us to relate these diagnostics to important physical parameters such as cloud mass or cluster age. Finally, we explore the implications of our models regarding the reliability of simple metallicity diagnostics, the properties of long-lived embedded clusters, and the role played by winds and supernovae in regulating \hii region and PDR line emission.


INTRODUCTION
Massive stars play a central role in the evolution of the interstellar medium (ISM) and the self-regulation of star formation. Feedback from these stars in the form of strong stellar winds, ionizing and non-ionizing radiation, and supernova explosions injects a large quantity of energy and momentum into the ISM and is one of the main processes that regulates the star formation rate on scales comparable to individual dense molecular clouds (see e.g. the reviews by Zinnecker & Yorke 2007, Krumholz et al. 2014and Klessen & Glover 2016. In addition, in all but the very nearest star-forming systems, observations of "star formation" are actually observations of young massive stars, either directly (e.g. in the far ultraviolet continuum) or more commonly indirectly, in the form of reprocessed radiation from gas (e.g. Hα, [C ii]) or dust (e.g. 24µm or far-infrared emission). Observations of line emission from the H ii regions and photodissociation regions (PDRs) created by radiation from massive stars are also widely used as diagnostics to estimate gas properties such as density or metallicity, particularly in extragalactic systems (see e.g. Kewley et al. 2001;Kauffmann et al. 2003;Kewley et al. 2013;Cormier et al. 2015;Sánchez et al. 2015).
Our understanding of stellar feedback, the ISM and how the latter responds to the former is based in large part on observations of nebular and molecular lines as well as dust emission powered by massive stars. These provide information about how radiation propagates in star-forming regions, the distribution of molecular gas, and the local heating of dust by starlight. Line intensities place constraints on gas density and emissivity, while line profiles inform us about line-of-sight gas motions, including turbulence, and large scale motion. Correctly interpreting the information provided by these various types of observations is of central importance to our understanding of star formation.
There is therefore a great need for models capable of predicting the properties of a wide range of observational tracers as a function of the gas distribution within a starforming cloud and the stellar population forming inside it. Unfortunately, this is a very difficult problem. One would ideally use high-resolution 3D simulations of massive star formation and feedback within molecular clouds as a basis for making these observational predictions. However, models of this type that include all of the relevant physics (photoionization, photo-dissociation, radiation pressure, winds and outflows, as well as supernovae) remain thin on the ground, with most 3D simulations only including some subset of these processes (see e.g. Dale et al. 2014;Rosen et al. 2016;Peters et al. 2017a;Haid et al. 2018, for some recent examples). Moreover, these models are computationally expensive, often limiting them to relatively low resolution, and making it completely impractical to use them to explore a large parameter space.
For this reason, the interpretation of Hii region and PDR diagnostics has often been carried out using models that do not attempt to self-consistently follow the growth of an Hii region and PDR around an assembling cluster of massive stars, but instead simply consider a broad grid of values for the gas density, radiation field strength and/or ionization parameter 1 , the composite spectrum of the stellar cluster, and often also the metallicity (see e.g. Kaufman et al. 1999;Kewley et al. 2001;Kewley & Dopita 2002, for some influential examples of this type of model). More recently the calculations have also systematically included variable cluster ages and time-varying spectral energy distributions (SEDs) (see Byler et al. 2017), while not solving for the coupling between cloud and cluster age. This approach is convenient and computationally efficient, but relies on the often unstated assumption that these parameters are not strongly correlated, such that it is meaningful to vary them independently of each other. As we will see later, this is an over-simplification that does not correctly describe the behaviour of real Hii regions.
One-dimensional dynamical models offer a useful compromise between computationally expensive 3D models and computationally efficient but physically unrealistic static models. Although these models have a long history (see e.g. Castor et al. 1975;Weaver et al. 1977), most make simplifying assumptions such as the neglect of gravity and/or radiation pressure that limit their applicability. In Rahner et al. (2017, hereafter Paper I), we introduced warpfield, a new 1D code for simulating the time evolution and structure of the stellar wind bubble, Hii region and PDR surrounding a cluster of massive stars. This code accounts self-consistently for the physics of stellar winds, supernovae, radiation pressure, ionization and gravity, allowing it to be applied to clouds of finite size that are not necessarily immediately disrupted by stellar feedback. It solves explicitly for the density structure adopted by the gas in response to the action of the feedback, with the help of the reasonable assumptions that the internal pressure of the feedback-blown bubble is much larger than the external pressure and that the shell surrounding the bubble is in quasi-hydrostatic equilibrium (Abel et al. 2005;Pellegrini et al. 2007Pellegrini et al. , 2011. Finally, the time-dependent approach adopted in warpfield also allows one to account for the evolution of the luminosity and composite spectrum of the stellar cluster as it ages. Although warpfield has proved to be a useful and flexible code for analyzing the evolution of feedback-affected regions around star clusters (see e.g. Rahner et al. 2018;Rugel et al. 2019, for some initial applications), it does not by itself allow us to predict the properties of observational tracers of these regions. In the present paper, we address this by coupling warpfield with state-of-the-art treatments of the microphysics of gas in Hii regions and PDRs (cloudy; see Ferland et al. 2017), and of line and continuum radiative transfer (RT) (polaris; see Reissl et al. 2016). This coupled model, which we refer to as the warpfield emission predictor, or warpfield-emp for short, allows us to make detailed predictions for the time-dependent line and continuum emission arising from the Hii region and PDR surrounding an evolving star cluster (potentially including on-going star formation). It accounts for both nebular and dust emission and also for the important effects of dust absorption on the emergent spectrum.
The structure of our paper is as follows. In Section 2, we review the basic features of warpfield and describe several recent updates we have made to the code. We then discuss how we couple warpfield to cloudy and polaris to construct the combined warpfield-emp framework. We also discuss how we model the SED from the central stellar cluster. In Section 3, we present an important application of this framework, an exploration of how the location of a star-forming region in a line diagnostic diagram such as the BPT diagram (Baldwin et al. 1981) varies as the region evolves. Section 4 compares warpfield-emp output with real observational data and also discusses some of the complications involved in doing so in a meaningful fashion. Section 5 examines whether we can use the location of our model star-forming regions in the BPT diagram to draw inferences about the initial and current properties of the host cloud and stellar cluster. In Sections 6 and 7 we discuss some important implications of our results regarding the applicability of simple metal abundance estimators and the nature of embedded clusters. In Section 8, we explain why our dynamical approach to modelling line emission from starforming regions is a significant improvement over previous approaches. Finally, in Section 9 we present a summary of our main results.

WARPFIELD-EMP
The warpfield-emp framework consists of three components: warpfield (Rahner et al. 2017), used to model the evolution and physical structure of the Hii region and PDR; cloudy (Ferland et al. 2017), used to determine the chemical state of the gas, the dust temperature, and the emission properties of the gas and dust; and polaris ), used to account for the effects of reddening on the emergent radiation. As all three codes are described in detail elsewhere, here we will only give a brief outline of each, with a focus on any modifications we have made compared to the previously published versions. We will also discuss how we couple the three codes together in the warpfield-emp framework.

warpfield
warpfield solves the equations of motion for a 1D spherical shell 2 evolving under the influence of stellar winds, supernova energy injection, radiation pressure and gravity. It assumes that the evolution of the region affected by feedback can be separated into several distinct phases, illustrated schematically in Figure 1. In Phase I, the shell is filled with hot gas from shocked stellar winds, and its expansion is driven by the pressure difference between the hot interior gas and the much colder ambient medium. During this phase, the effects of gravity and radiation pressure are negligible, leaving one a relatively simple differential equation to solve for the shell motion (Weaver et al. 1977;Bisnovatyi-Kogan & Silich 1995;Rahner et al. 2017). Phase I comes to an end once the bubble loses its hot gas, either because the gas cools, which occurs after a cooling time t cool , during shell-fragmentation , or because the bubble "bursts", allowing the hot gas to escape through lower density channels (Rogers & Pittard 2013). In our 1D treatment, we cannot solve explicitly for effects like the bursting of the bubble, which are influenced by the 3D structure of the cloud, and so by default we assume that the bubble bursts once the shell has swept up the entirety of the initial cloud, i.e. once the shell radius is equal to 2 Although our 1D approach is primarily motivated by the desire to avoid the computational costs of a full 3D approach, which are too high to allow one to properly survey the relevant parameter space, there is also a physical justification underpinning this. The internal structure of H ii regions traced by line ratios reveals them to be largely regular, mostly circular objects when studied in ionization parameter mapping line ratios (Pellegrini et al. 2012) in the Magellanic clouds or NGC 628 (Rousseau-Nepton et al. 2019). Thus despite the variability in initial conditions, the cumulative effects of all feedback seem to wash out inhomogeneities and often produces spherical objects that should be described well by our 1D models. the initial cloud radius. This occurs at a time tsweep, which we can readily determine from our 1D treatment. Phase I therefore comes to an end after a time t = min(t cool , tsweep). 3 Following Phase I, the shell either proceeds to Phase II (if t < tsweep) or directly to Phase III (if t tsweep). In Phase II, the shell continues to sweep up material from the surrounding molecular cloud. The expansion of the shell is driven by the ram pressure exerted on the shell by stellar winds and supernovae, and by the radiation pressure acting on the shell, but at this point we also start to account for the counteracting effect of gravity, both due to the gravitational attraction of the central cluster on the shell and also the shell's own self-gravity. Phase III is very similar to Phase II in terms of the force balance acting on the shell, with the key difference being that during Phase III, the original cloud has been entirely swept up, with the shell now assumed to be expanding into the low density warm neutral medium (WNM).
There are two possible fates for the shell. If gravity begins to dominate, which often happens as the central stellar cluster ages and feedback becomes less effective, then the shell will begin to recollapse. In the original warpfield implementation (Rahner et al. 2017), we simply terminated the calculation when the shell had recollapsed to a radius of 1 pc. In warpfield-emp, however, we adopt instead the approach used in Rahner et al. (2018) and assume that recollapse triggers a new burst of star formation, which for simplicity is assumed to have the same star formation efficiency as the original burst. The other possible fate is continued expansion without recollapse, which will occur once the kinetic energy of the shell exceeds its gravitational potential energy. In this case, we follow the expansion of the shell until its maximum density falls below 1 cm −3 for an extended period of time (more than 1 Myr). Following this, we con-sider the shell to have dissolved. In the analysis in Section 3. we also ignore models where the shell has expanded to more than 200 pc since at this point the influence of environmental effects such as galactic shear or interactions with other clouds would start to become important -issues which are currently ignored in warpfield.
As part of warpfield, we solve not only for the dynamics of the shell, but also for its structure, and the density gradients driven by radiation pressure, as well as the natal cloud structure. This strongly influences the effectiveness of radiation pressure driving. To do this, we assume that the ionized and neutral/molecular portions of the shell are in quasi-hydrostatic equilibrium, with the equation of state given by Abel et al. (2005) and Pellegrini et al. (2007). By default, warpfield accounts only for thermal pressure support when solving for the shell structure. However, the code offers the option of also including the effects of magnetic and/or turbulent pressure support, as discussed at greater length in Appendix B. In the results presented in this paper, we do not make use of these options in order to minimize the number of free parameters.
The four main input parameters to any warpfield model are: the mass of the gas cloud (M cl ); the initial mass of the stellar cluster, (M * ), which can either be specified directly or in terms of a star formation efficiency , in which case M * = M cl ; the density of the cloud (n cl ); and the metallicity. In the models presented later in this paper, we keep the metallicity fixed at the solar value, but explore the effects of varying M cl , and n cl . We assume that the stars within the star cluster follow the Kroupa initial mass function (Kroupa 2002) where the number of stars N with initial mass M is given by (1)

Cloudy
Cloudy is a non-local thermodynamic equilibrium (NLTE) spectral synthesis and plasma simulation code developed over the past three decades by G. Ferland and collaborators. In warpfield-emp, we use the latest major version of Cloudy, version C17, described in detail in Ferland et al. (2017). We use Cloudy to determine the emissivities of a large set of emission lines as a function of position within the shell and the surrounding undisturbed cloud. Cloudy also provides us with a description of how the corresponding opacities vary within the shell and the cloud. This data is then passed on to polaris to produce the final emission line maps, as described in Section 2.3 below. For shells in Phase III, the procedure is straightforward. In this case, the shell has swept up the entirety of the surrounding cloud and so we need only consider emission from the shell itself. To do this, we pass the following information to Cloudy: the frequency-dependent flux from the stellar cluster that is incident on the inner edge of the shell, the radial distance from the central cluster to the shell, the total pressure (ram plus thermal) acting on the inner shell boundary, the velocity of the shell, the total shell mass, and the metallicity. With the exception of the frequency-dependent flux, which we generate using the time-dependent cluster model discussed in Section 2.4 below, all of these parame-ters either come directly from our warpfield model for the shell or from the initial conditions for the calculation.
Given this information, Cloudy then solves for the structure of the shell by breaking it up into a series of thin concentric shells ("zones" in the terminology used in the Cloudy documentation). The thickness of each of these zones is determined internally by Cloudy and is set by the requirement that the physical conditions across each zone must be close to constant. Consequently, the zone thickness can vary significantly as a function of distance from the cluster -for example, the zones in which the ionization front are located must be very thin in order to properly capture the transition from ionized to neutral gas. The density in each zone is determined by Cloudy, using the same assumption of quasi-hydrostatic equilibrium as in warpfield. 4 cloudy also solves self-consistently for the thermal structure of the gas, its chemical composition, and the emissivities and opacities for a large set of different emission lines.
For shells in Phase I or Phase II, which have not yet swept up all of the natal cloud, there is one additional step to the procedure. Once we have determined the shell structure and emission, we then carry out a second Cloudy calculation to determine the emission from the natal cloud. In this case, the incident flux is taken from the output of the Cloudy calculation for the shell. 5 The cloud is assumed to have zero bulk velocity relative to the cluster, and to have a turbulent velocity dispersion set by the requirement that the cloud be initially in virial equilibrium. The other parameters that are required -the mass, density and size of the cloud, the location of the outer edge of the shell, and the metallicity -come directly from the warpfield calculation or the initial conditions. As in the case of the shell calculation, Cloudy breaks the cloud up into a series of zones and computes the emissivities and opacities for each zone.

Polaris
The time-dependent evolution of density and column density is large in our models. The resulting internal extinction by dust plays a major role in determining the emergent lines intensity and therefore the ratios. For example, the ratio of the ionizing photon production rate, Q0, to the emergent Hα flux varies over 9 orders of magnitude between ultracompact H ii regions with AV 20 and unembedded H ii regions with AV 1. Our model predicts regimes where the illuminated portion of the shell that is the H ii region is heavily obscured by surrounding material, meaning that it is possibly unobservable in many tracers. Further complicating the picture, radiation making it out of the H ii region into the lower density ambient cloud can produce low ionization emission that can dominate over the emergent light from the H ii region (see Figure 2). It is therefore extremely important to account for the impact of dust extinction. To do this, we Figure 2. A schematic of an embedded H ii region, surrounded by a partially ionized cloud. The net emission in Hβ, more blue and attenuated by dust than Hα, may be dominated by the static component, while Hα emission from a H ii region contributes measurably to the net flux from the cloud. In such a scenario, common to many of our models, the Balmer decrement is affected by differential attenuation and therefore, it has no simple physical meaning and any reddening correction derived from it may be highly inaccurate. make use of a modified version of the polaris RT code 6 .
polaris was originally developed to model polarized dust continuum emission from 3D dust distributions using the Monte Carlo method (see e.g. Reissl et al. 2017Reissl et al. , 2018 for some example applications), but more relevant for our purposes is its ray-tracing mode. It has been used for various applications, such as state-of-the-art treatments of dust grain alignment physics Reissl et al. 2016Reissl et al. , 2017Seifried et al. 2019;Reissl et al. 2018), atomic and molecular line RT including the Zeeman effect (Brauer et al. 2017;Reissl et al. 2018), as well as synchrotron polarization and Faraday rotation (Reissl et al. prep) on variable 3D grid geometries (adaptive octree, spherical, cylindrical, and native Voronoi).
In ray-tracing mode, polaris solves the multi-frequency equation of RT on a set of parallel rays passing through a 3D grid, self-consistently accounting for the absorption and emission of radiation by dust. When used for line transfer, it can either compute the line emissivities internally, using e.g. the local thermodynamic equilibrium (LTE) or large velocity gradient (LVG) approximations to determine the local levels population, or it can adopt pre-computed values per grid cell supplied by the user. In warpfield-emp, we use the latter option, supplying the emissivities computed by cloudy as an input to polaris. In addition, we also supply polaris with the line and continuum opacities determined by the cloudy calculation 7 . 6 http://www1.astrophysik.uni-kiel.de/~polaris/ 7 In practice, for the diagnostics examined in this paper, the line polaris supports several different types of 3D grid, but in warpfield-emp the obvious choice is a 3D spherical grid, as in this case the transfer of information from the 1D model to the 3D grid is trivial. When setting up our spherical grid, we ensure that the length of each radial bin dr matches the thickness of the corresponding Cloudy zone, allowing us to map data directly from the zones to the radial bins. This required us to modify polaris to support spherical grids with arbitrarily-sized radial bins and potentially with large changes in thickness between bins. This is important for resolving the narrow emission line region near the ionization front, without requiring us to use a prohibitively large number of radial bins. We then carry out a ray trace along a set of parallel rays passing through the grid, allowing us to project the 3D model onto a 2D detector. When projecting from 3D onto our 2D detector, we use a uniform pixel matrix. Here, we make use of polaris's recursive sub-pixelling technique, which splits up rays as necessary to ensure that every cell in the 3D grid is sampled by at least one ray, thereby ensuring that we do not lose information from regions smaller than the size of the pixels. Finally, for carrying out the ray-tracing, cubic spline interpolation is used to interpolate values from the grid onto the rays, and the Runge-Kutte-Felhberg (RKF45) solver with an inbuilt error and step size correction is used to solve the RT equation along the rays. Since the radial distances dr between adjacent grid cells are already optimized by cloudy we limit the RKF45 solver to refine its step size only down to 0.1 × dr without losing accuracy. This enables us to ray-trace the enormous number of cloudy models within a reasonable time with an error below 1 %.
With this approach, we can produce maps of velocityintegrated surface brightness of any line or continuum tracers, provided that they are included in the cloudy model. Alternatively, we can produce full position-position-velocity (PPV) cubes of the emission line. In the latter case, the line profiles are calculated using a specified grid of velocity channels and accounting for the line-of-sight projected velocity, turbulence and thermal broadening of each line. In calculations that include both the shell and the remaining natal cloud, we make the simplifying assumption that there is no line overlap between the two components when the shell expansion velocity is much greater than the turbulent velocity dispersion of the cloud (vexp σ turb ). Optionally, the ray tracing can be perform with no dust opacity, recovering what would be the line emission without internal extinction.
In Figure 3, we show an example of the kind of output that we can produce with warpfield-emp. Each quadrant shows a section of the map we obtain by taking the ratio of the surface brightness of some observational tracers. The ratio [O iii]/Hβ, shown in the lower quadrant, traces regions with a large flux of energetic ionizing photons and hence is high throughout much of the projected bubble. However, it drops off dramatically toward the edge of the bubble as the ionizing flux drops, and the internal density increases. The opacities are very small and the dominant contribution comes from the dust. However, our method does not require this to be the case and can also be applied to produce maps of tracers with high line opacities, such as CO. ratio [S ii]/Hα, shown in the top quadrant, has very different behavior. It is low within most of the low density bubble (as most of the sulphur present has been ionized to S 2+ ), but increases dramatically as one nears the dense shell. This behavior is characteristic of optically thick H ii regions, and results in a pronounced increase in the [S ii]/[O iii] ratio at the H ii region boundary that is observed in many resolved H ii regions (see e.g. Pellegrini et al. 2012). The layer with enhanced [S ii]/Hα is thicker than the shell itself, which is only ∼ 0.05 pc wide at this point in the evolution of the bubble. We also see similar morphology in the predicted [N ii]/Hα flux (right) and [O ii]/Hβ flux (left). The region with elevated [N ii]/Hα is thicker than that with elevated [O ii]/Hβ, due to the higher ionization potential, even though the total intensity of [O ii] is higher owing to the higher abundance of oxygen relative to nitrogen. It is important to note that much of the details of radiation transfer, such as the inner boundary, are difficult to observe directly.

Individual Cluster Evolution
In order to allow us to focus in this paper on the effect of a time-dependent, dynamic cloud structure on predictions for emission line diagnostics, we adopt a very simple star cluster model for the simulations presented in this paper.
We assume that a given burst of star formation forms a star cluster with mass M cl , where is the star formation efficiency and M cl is the mass of gas in the cloud. The star formation efficiency and the initial mass of the cloud are input parameters of our model. For bursts of star formation occurring following re-collapse of the shell, we adopt the same star formation efficiency as for the initial burst, and compute the mass of gas in the cloud by subtracting the mass of stars already formed from the initial cloud mass. The stellar population of the cluster formed by a given burst of star formation is then simply scaled from a fully sampled M = 10 6 M cluster. As in our previous warpfield papers, we use the time-dependent evolution of a cluster calculated from Starburst99 8 (Leitherer et al. 1999(Leitherer et al. , 2014, using the Geneva stellar evolution tracks for rotating stars (Ekström et al. 2012). In a subsequent paper, we will introduce the complexity of stochastically sampled clusters when we consider populations of star-forming regions, but for now we focus only on predictions and degeneracies introduced with a fixed cluster template.

On-going Star Formation
While an individual stellar cluster may be considered "simple", the ionizing spectral energy distribution (SED) of our models is not. Unlike typical photoionization grids employing a single assumed cluster SED, our SED not only evolves in time with our physical conditions, but the inclusion of subsequent bursts of star formation following re-collapse of the shell leads to complex stellar populations, composed of multiple massive clusters with intermediate ages (few Myr), dominated at times by O-and/or B-stars. Normally the clusters with O-stars dominate the ionizing energy regime, but that is not guaranteed. Multiple older clusters with ages ranging from 2-10 Myr can be the dominant source of radiation in other frequency bins. This is especially true below the hydrogen ionization limit, and can greatly affect the ratios of low-to-mid ionization potential lines to high-ionization potential emission.
As an example, we show in Figure 4 the case of a cluster that has undergone three separate bursts of star formation: an initial burst, plus two subsequent bursts associated with shell re-collapse. At short wavelengths -particularly shortwards of 912Å -the SED of the cluster is dominated by the contribution from the youngest stellar population, while at long wavelengths, the two older stellar populations dominate.

TIME-DEPENDENT DIAGNOSTIC DIAGRAMS
One of the major strengths of the warpfield-emp framework is the ability it gives us to examine the time-dependent behavior of the emission from an evolving H ii region and PDR in a physically self-consistent fashion. As a star cluster ages, its SED changes (Leitherer et al. 1999), and so the relative strengths of the emission lines associated with the ionized gas also change. However, the strengths of these lines are also sensitive to the density of the emitting gas, which changes as the H ii region evolves, as well as other variables such as the distance of the ionized gas from the cluster. Moreover, these variables are correlated, since the shell density, radius etc. depend to a great extent on the time history of feedback from the central cluster. Therefore, approaches in which these variables are varied independently are likely to yield misleading results.
As an example of the power of the warpfield-emp approach, we examine the example of a pair of widely-used observational diagnostics of ionized gas, the ratios of the One of the most common uses of these line ratios is to identify the mechanism responsible for ionizing the gas. Baldwin et al. (1981)  To explore how we expect the [O iii]/Hβ and [N ii]/Hα ratios to change as a function of time in evolving H ii regions, we have carried out an extensive series of simulations using warpfield-emp. As summarized in Table 1, we have examined models with three different initial densities, six different star formation efficiencies, and cloud masses ranging from M cl = 10 5 M to 10 7.25 M in 0.25 dex increments, yielding a total of 180 different models. For simplicity, in this initial exploration of parameter space, we do not vary the metallicity, keeping it fixed at Z = Z in every run. This helps to limit the number of models we need to run, and also allows us to explore how much variation we expect in the [O iii]/Hβ and [N ii]/Hα ratios at fixed metallicity.
To limit the model output data volume and avoid unnecessary computational work, we do not re-run cloudy and polaris after every warpfield timestep. Instead, we delay re-running cloudy until the physical conditions in the cloud or the properties of the star cluster have changed sufficiently to result in an appreciable change in the emission. In practice, this means that we re-run cloudy whenever the shell density at the inner shell boundary, shell radius, shell mass or the ionizing photon flux change by the fractional amounts shown in Table 2. We also recalculate models when the warpfield model changes evolutionary phases, or if more than 0.5 Myr has elapsed since the last calculation.
In Figure 5, we show the results from an example model with M cl = 10 5.00 M , nH = 100 cm −3 and SFE = 1%. In grey scale we show the full distribution of our models. Further examples are shown in Appendix C for different cloud properties. Both observables, such as the Hα luminosity, and intensive properties, such as the density at the bubble/H ii region boundary, evolve rapidly with time. At early times the cluster luminosity is relatively constant, and this evolution is driven by the change in internal extinction. The density plotted is the density at the bubble boundary, and is set by the internal hot wind pressure, or after the bubble has burst, cooled, or fragmented, the direct momentum of the free streaming wind. For the ionization front density, as with all radiation-pressurized nebulae, there is a relationship between density diagnostics and the photon flux incident on the cloud. This produces a fundamental relationship between star-forming events and density, but one that is less sensitive to initial conditions and more a reflection of the current nebular structure and cluster age. Thus, different trends with ionized gas densities identified in nearby galaxies in Herrera-Camus et al. (2016) should be studied in the context of full population synthesis, and are beyond the scope of this immediate paper. It is sufficient to say that these trends are related to, but not caused by environment.
In the second row, we plot the warpfield evolutionary phase and ages of the models. We see immediately that both line ratios change substantially as the cluster and the shell evolve, particularly the [O iii]/Hβ ratio. To help understand this we introduce a metric to quantify the relationship between geometry and emission, r obs,fraction , to represent the radius from the bubble boundary where the cumulative luminosity of a particular observable reaches a certain fraction of the total the model produces.
Two uses of the ratio are shown in the last two panels of Figure 5, where we show the values of rHα,90 − rHα,10 which is a proxy for the observable thickness of the H ii region, as traced by Hα, a better way of quantifying whether the shell is "thick" (in which case the value of the thickness is large compared to the inner radius of the bubble) or "thin" (in which case the thickness is small compared to the size of the bubble). We can also use [O iii] and Hα because their strong differences in ionization potential guarantees that we are tracing the extremes of the H ii ionization structure. . Example of a cluster with a complex star formation history. In this model, which has a cloud mass M cl = 10 6 M , a cloud density n H = 500 cm −3 and a star formation efficiency = 0.02, the cluster at the time depicted has undergone three bursts of star formation -an initial burst, and two subsequent bursts associated with shell recollapse. In the top panel, we show the spectral energy distributions (SEDs) of the three distinct stellar populations, at a time when their ages are t = 0.4, 6.3, and 13 Myr, respectively. The SED of the oldest cluster at time t = 0 is shown with the gray line. The total mass of the cluster (in log units) and the masses of the individual stellar populations are listed in the inset. We see that the ionizing luminosity and the B-band are dominated by emission from the youngest cluster, while the light at red wavelengths and in the IR is dominated by the contribution of the older clusters. This is also highlighted in the bottom panel, where we show the fractional contribution to the total luminosity at each wavelength made by each separate stellar population.
At early times, during the first 2 Myr, the shell expands at high velocities, v > 10 km s −1 , gradually slowing down as it sweeps up more and more material. In the particular example shown, the fast growth of the shell coincides with an increase in [O iii]/Hβ at very early times. The increase comes as the adiabatic expansion decreases the internal pressure, leading to a decrease in the density of the ionized gas in the shell. This density is initially very high because the high internal pressure in the bubble strongly compresses the surrounding ionized shell. As the shell accelerates, the internal volume increases, leading to a decrease in the density of the shell and an increase in the [O iii]/Hβ ratio.
At later times, the emission rate of ionizing photons starts to decrease when the most massive stars die, and the [O iii]/Hβ ratio drops again. The shell expansion velocity, which has decreased to ∼ 5 km s −1 , increases again as the ram-pressure of the first SN explosions drives the gas away from the star cluster. Shortly afterwards, the nebular emission reaches its maximum Hα luminosity of ∼ 10 38 erg s −1 . After approximately 10 Myr, as the shell extends over an ever increasing volume, the density in the shell drops so low that it becomes easier to produce [O iii] again even though Q0 continues to drop. However, we note that the evolution of the presented example model should not be taken to be representative of all star-forming regions. Objects with a stronger starburst, where the cloud is quickly destroyed by feedback, and objects with a lower SFE where consecutive star formation events will occur, can -and often will -have a very different evolution in a BPT diagram, as the additional examples shown in Appendix C make clear.
A final note on recollapsing objects: We often find that the density of the shell in a recollapsing bubble decreases even as the radius becomes smaller. This is a reflection that the feedback is decreasing faster than the gas can compress, manifesting as a thicker H ii region with higher filling factor and lower average density. In some time series of more massive clouds, this balance can lead to an increase in [O iii]/Hβ.
Although informative, plots such as those in Figure 5 that show where the emission from a single star-forming cloud falls in the BPT diagram throughout the lifetime of the cloud are also potentially misleading, since not all evolutionary stages are observable. Early on, before the whole of the natal cloud is swept up into the shell, the foreground dust extinction due to the undisturbed portion of the cloud can be large. This has little impact on the location of the points in the plot, since the values of the [O iii]/Hβ and [N ii]/Hα line ratios are insensitive to extinction. However, the absolute brightness of the lines is strongly affected and suffi- The line shows the evolutionary track, while the circles represent the predictions at each non-uniform time step. In each panel color represents an observable, a cloud property, or property of the system. For each star formation event the first data point after the star formation event (at 0.1 Myr) and the occurrence of first supernova (at 3.6 Myr after cluster formation) are marked by a star symbol and cross, respectively, with a number indicating which sub-cluster the event belongs to. The end-point of the model is indicated by an arrow. Symbols outlined in red indicate times at which the system is unlikely to be observable, given typical sensitivity limits for extragalactic observations. For reference, we also indicate the criteria proposed by Kewley et al. (2001) (dotted line) and Kauffmann et al. (2003) (dashed line) for distinguishing between star-forming regions and AGN. The results from our full set of models at every output time are shown in the gray-scaled hexbin histogram in the background. Contributions to this histogram are weighted by the size of the time-step. From left to right and top to bottom, the panels show: Emergent, dust-attenuated Hα luminosity; visual extinction (A V ) along the line of sight, computed assuming R V = 3.1, total H nuclei number density at the wind-bubble-H ii region interface; warpfield evolutionary phase, with 0 corresponding to re-collapse; age of system from the beginning of the simulation; age of the youngest star cluster (accounting for re-collapse); radius of the bubble/shell boundary in pc; difference of the radii (as measured from the bubble/shell boundary) containing 90% and 10% of the integrated Hα luminosity (r Hα,90 , r Hα,10 ); and lastly the shell velocity, with positive values corresponding to expansion and negative to re-collapse. cient dust may render them undetectable given any plausible amount of observing time. At late times, on the other hand, the H ii region becomes faint not because of dust extinction but simply because by this time, most of the massive stars powering the emission have already died. In addition, the emission that is present is spread over an increasingly large area of sky, owing to the expansion of the shell, further reducing its surface brightness. In order to compare results from warpfield-emp with observations of real H ii regions in a meaningful fashion, it is necessary to account for both of these effects. We discuss how we do this in Section 4 below.

COMPARISON WITH OBSERVATIONS
There are two different ways in which we could consider comparing the results of our warpfield-emp models with real observational data. The first involves the comparison of models and observations of a single region. For example, in previous work we have found warpfield models that are a good match to the current evolutionary state of the H ii regions observed in 30 Doradus in the Large Magellanic Cloud ) and the W49 star-forming region in the Milky Way (Rugel et al. 2019). Using warpfield-emp, it is straightforward to make predictions for the line and continuum emission from these regions that can be compared to the observed values. In this case, little post-processing is required to convert the surface brightness maps produced by warpfield-emp into a form that can be directly compared with the observations. The maps must be smoothed to match the angular resolution of the observations and noise must be introduced at a level that is consistent with the noise in the real data. As the appropriate values to choose here will be specific to the particular observations of interest, warpfield-emp itself does not account for either of these post-processing steps, which must therefore be carried out by the user. This is beyond the scope of the present paper.
The other type of observational comparison is one between a population of objects and a range of different warpfield-emp models. This is what we focus on here. It is important, as it tells us how well our models succeed in reproducing the diversity of different H ii regions and winddriven shells observed in real galaxies. A proper comparison along these lines requires a population synthesis approach to reproduce the observed density of models at each value of [N ii]/Hα and [O iii]/Hβ ratios, since we cannot assume that we are observing all of the different star-forming regions at the same time in their history. Doing this fully selfconsistently lies beyond the scope of the current paper, although we will return to this point in future work. However, we can approximate the observed distribution in a diagnostic such as the BPT diagram by appropriate selection and weighting of our model predictions, as we describe in more detail below.
The set of observational data that we use for the purposes of comparison is a catalog of H ii regions based on tunable filter observations of NGC 628 from Rousseau-Nepton et al. (2018). That catalog provides line fluxes integrated over the entire extent of each H ii region. We are using a version of the catalog where no reddening correction has been applied, which is significant for comparisons involving [O iii]. Direct comparison of our model predictions with a catalog such as this requires some care, as the survey characteristics and catalog classification methods affect the number and properties of the regions included in the catalog.
To take a simple example, any optically-selected catalog of H ii regions has a clear bias against including young, highly embedded objects, since these will have high extinctions and hence will often have optical line fluxes below the detection limit of the survey. Failing to account properly for effects such as this can therefore lead to apparent discrepancies between models and observations that are purely a consequence of the observational selection effects.
In order to make a meaningful comparison between our set of models and the the Rousseau-Nepton et al. (2018) H ii region catalog we employ the following procedure, outlined in Figure 6, to produce Figure 7. Note that in what follows, we distinguish between a warpfield-emp "model" (the time history of the emission produced for a single combination of cloud mass, SFE, cloud density etc.) and a "model output" (the results of a warpfield-emp model at a particular output time).

• Step 1:
We produce a list of the H ii regions in the Rousseau-Nepton et al. (2018) catalog ordered by increasing observed Hα luminosity, L(Hα). We then divide up the H ii regions from the list into three bins, the first containing the first 25% of the regions in the list (i.e. the 25% of H ii regions with the lowest observed Hα luminosity), the second containing the next 50% and the last containing the final 25% (row 1, column 1 in Figure 6). In practice, this results in a set of bins with boundaries 35.8 < log 10 L(Hα) < 36.7, 36.7 < log 10 L(Hα) < 37.5 and 37.5 < log 10 L(Hα) < 39.4, with Hα luminosities in erg s −1 . Binning the H ii regions in this way allows us to (approximately) delineate regions with different masses and ages. Each bin corresponds to a different row in Figure 7. The H ii regions from the Rousseau-Nepton et al. (2018) catalog are reddened not only by dust within the region itself (which is accounted for in our warpfield-emp models) but also by foreground dust within NGC 628 and the Milky Way, which is not accounted for in our models. When selecting the set of model outputs to compare to each observational bin, we therefore allow for a range of foreground extinctions 0 < AHα,max < 1.5 mag (row 2, column 1 in Figure 6). Using this range, we identify all model outputs where the predicted emergent Hα luminosity falls within the range Lmin < L(Hα) model < Lmax × 10 A Hα,max /2.5 , where Lmin and Lmax are the minimum and maximum Hα luminosity included in the bin. • Step 2: We next divide the observational data into low, medium and high uncertainty subsets, using the size of the 1σ error in the . This gives us a description of how the error in the line ratios varies as a function of the ratios for the H ii regions in each bin. Note that this approach is more accurate than adopting a single representative error value for each ratio in each bin, as it accounts for the substantial covariance between the errors and the line ratios.
• However, if we want to compare these values in a meaningful way with the observed values, we need to account for the uncertainty in the observations. Essentially, what we would like to do is to assign each data-point a set of line ratio errors that are comparable to those in the observational data. This is made more complicated by the fact that the size of these errors differs in the different observational subsets and also varies systematically as a function of the line ratios within each subset. Therefore, rather than assigning each model output the same error values, we instead make use of the 2D spline fit described in Step 2. In each Hα luminosity and [O iii]/Hβ error bin, we run through the set of model outputs associated with that bin. For each model output, we determine whether its predicted line ratios lie within 3σ of at least one observational data-point. If they do, then we use the spline fit to compute the appropriate errors that we should assign to the model data-point. If they do not -i.e. if the model output lies in a part of the [O iii]/Hβ-[N ii]/Hα parameter space that is not represented in the observational data -then we cannot safely use the spline fit. In this case, we instead simply assume that the errors in these model outputs are equal to the median error of the observed points in that particular Hα luminosity and [O iii]/Hβ error bin. Finally, to produce the smooth distribution of simulated data illustrated by the color maps in Figure 7, we represent each model output data-point and its associated errors as a 2D Gaussian and carry out a weighted sum of these Gaussians. 10 9 We use [O iii]/Hβ rather than [N ii]/Hα to distinguish the different subsets as the observational error in the former is generally higher than that in the latter. 10 In practice, we do this by Monte Carlo sampling the individual As weighting factors, we take the timestep corresponding to that model output 11 in units of 10 5 yr, as well as a factor accounting for the probability of finding a cluster of the particular mass associated with that model output. This latter factor is computed assuming a cluster mass function with a slope of α = −2.0. Note that by construction our set of model outputs only includes results for shells with radii r 200 pc, since we assume that on larger scales the effects of galactic shear will tend to disrupt the shell, and hence terminate our warpfield calculations for shells that reach this size. We therefore confirm that all of the H ii regions in the utilized observational dataset have sizes smaller than this. In practice, this restriction turns out to be unimportant, as almost all H ii regions observed in NGC 628 have radii r 200 pc.
Overall there is good agreement between the models and the observations. Most of the observed H ii regions lie on a locus in the [O iii]/Hβ-[N ii]/Hα plane that is also densely populated by model results, meaning that many of our models spend a large fraction of the time for which they are observable in this region in the plane. Moreover, comparison with the high signal-to-noise data-points shows that the models not only fall in the right location in the plane, they also reproduce the 0.25 dex scatter seen in the [N ii]/Hα ratio at constant [O iii]/Hβ. Note that for this subset of observational data, the scatter is not simply explainable as a consequence of the statistical errors in [N ii]/Hα, as these are small. Interestingly, we reproduce this scatter despite our models having a constant N abundance, suggesting that it is primarily driven by variations in the temperature and excitation of the [N ii]-emitting gas, rather than by variations in the metallicity.
Nevertheless, we also see that there are regions in the plane where we have model results that correspond to few or no observed H ii regions. For example, at log [N ii]/Hα ∼ 1 and high log [O iii]/Hβ ( 0), we find a clear concentration of model outputs but only a few actual H ii regions. As we discuss in Section 5, models producing outputs in this part of the parameter space typically have high initial cloud densities, so the lack of observed H ii regions in this area of the plot may simply be telling us that there are few giant molecular clouds (GMCs) with such high initial densities in NGC 628. This would be consistent with the relatively small cloud-scale molecular gas densities that one infers from ALMA observations of GMCs in NGC 628 (Sun et al. 2018), although one should bear in mind that the latter values are likely biased low by the effects of beam dilution.
The second main region of interest are low ionization objects, with log[O iii]/Hβ −1.0 and log[N ii]/Hα ∼ −0.6. A few of our models do in fact populate this part of the observational plane. However, it remains under-sampled as it cor-responds to calculations with M cl 10 5 M and with star cluster masses below a few 1000 M . For such small clusters, our assumption of a fully-sampled IMF is no longer valid, and so reproducing the observational data in this region will likely require us to consider models with stochasticallysampled IMFs, a topic which is outside of the scope of this paper.

Dependence on Initial Cloud Properties
One interesting application of warpfield-emp is to explore whether the location of star clusters in the BPT diagram can tell us anything directly about their properties, or the properties of the cloud in which they formed. To this end, in Fig. 8 we show where models with different values for the main input parameters (initial cloud density, cloud mass, and input SFE) fall in the BPT diagram. In each hexbin in the plots we present the mode, i.e. the most frequent value inside that bin. This should not be confused with the actual observational probability, since our grid does not account for the fact that some input parameters -e.g. lower cloud masses -occur more often than others. We do, however, account for the use of adaptive time steps by weighting each point in the BPT-diagram by the length of its corresponding time step. 12 As can be seen in Fig. 8, the initial cloud density is not a good indicator of where a model will fall, even after pruning models due to observability. There is a slight tendency for H ii regions born in an initially low density cloud to populate the top left corner of the BPT diagram, i.e. low [N ii]/Hα and high [O iii]/Hβ, while H ii regions born in high density clouds tend to be found in the low [N ii]/Hα and low [O iii]/Hβ tail of the distribution. However, we also see that there are many high density models that spend part of their time in the low [N ii]/Hα and high [O iii]/Hβ region of the plot, and many low density clouds in the low [N ii]/Hα and low [O iii]/Hβ tail. This behavior is a consequence of the time dependence of the emission line ratios that we already discussed in Section 3. Because of this time dependence, a given cloud model does not produce line ratios that always remain in the same part of the BPT diagram, but instead yields values that move around over time.
We find similar behavior if we vary the other two input parameters, cloud mass and SFE. Clouds with high SFEs tend to be found at the top of the diagram, while very massive clouds mainly fall around the inner bend, but in each region of the plot there is also considerable scatter. Therefore, if the only information we have available on a given cluster and cloud is its location in the BPT diagram, we cannot conclusively infer the cloud or cluster mass, or the mean density of the cloud. The best we can do is to make a probabilistic statement, e.g. that an H ii region located at the top of the plot is more likely to be located in a low mass cloud than in a high mass cloud. Converting this qualitative 12 Roughly speaking, the simulation time step corresponds to the time a model spends inside a bin. This is only true as long as model does not move more than one bin size on the diagnostic diagram during one time step. Figure 8. Top: BPT diagram for our set of H ii region models, colour-coded by the value of the initial cloud density most likely to result in emission at that point in the diagram. Specifically, in each bin we show the time-step weighted mode, i.e. the value of n 0 that produces the set of models that collectively spend the most time in that bin. Middle, bottom: the same, but for the cloud mass and input SFE, respectively. Note that in all three panels, objects are only included if they have Hα luminosities exceeding 10 35.8 erg s −1 and sizes smaller than 200 pc, as in Section 4. statement into a quantitative statement and accounting for the fact that in reality the line ratios for any given H ii region are not known precisely is an important task, but one which lies outside of the scope of this preliminary investigation (see Ardizzone et al. 2018 for a first assessment).

Dependence on Evolved Cloud Properties
As well as looking at the connection (or lack thereof) between the initial conditions for each model and its later location in the BPT diagram, we can also investigate how this position correlates with properties defining the instantaneous conditions within the H ii region, such as the age at which we observe it. In Figure 9, we show BPT diagrams similar to those in Figure 8, but color-coded (a) by the age of the youngest cluster; (b) the ionization parameter, U , measured at the inner shell radius; (c) the size and (d) the expansion velocity of the shell; (e) the net star formation efficiency SFEtot; and (f) the time-averaged star formation rate (SFR). The net star formation efficiency is defined here as the ratio of the total stellar mass formed within the cloud to the initial cloud mass. In each case, we show the median value in each hexbin. In the cases where the feedback-driven shell never recollapses, the net star formation efficiency plotted in panel (e) is the same as the star formation efficiency SF that we specify as an input parameter. However, in cases where the shell collapses one or more times during the lifetime of the cloud, we form additional stars in a burst after each re-collapse, and in this case the net star formation efficiency is larger than the input value. Finally, the timeaveraged star formation rate shown here is defined as where M cluster,i is the mass of each individual cluster, with a total number of clusters equal to one more than the number of re-collapses, and ∆t0 is taken as the total age of the model. Unlike SFRs derived from large ensembles of objects, where it makes statistical sense to average the SFR over a given time frame (often taken to be 10 Myr; Kennicutt 1998), the values we report here assume that each individual SF event is instantaneous. In contrast to Figure 8, several clear correlations are visible in Figure 9. There is an obvious relationship between the age of the youngest cluster, the value of the ionization parameter and the position of the cloud in the BPT diagram, which we discuss in more detail below. There is also a noticeable but somewhat weaker correlation between the size of the shell and the value of the This is easily understood as an evolutionary effect. The size of the feedback-driven shell is correlated to the age of the youngest cluster. Compact expanding shells are generally found around young clusters, while star clusters whose feedback has driven the shell to radii 100 pc tend to be old. Therefore, objects found in the lower half of the BPT diagram are those with large shell radii. In galaxies with low shear, where shells larger than 200 pc remain coherent, objects with very large shells would also inhabit the region around the Kauffmann bend but would be hard to observe due to their low surface brightness. Note also that although compact shells can be found around older clusters in cases where the shell is collapsing, these objects are generally not observable. This is clear from panel (d) in Figure 9, where we see that only a few of the H ii regions that pass our observability cuts are collapsing. The majority are expanding, with velocities ranging from a few up to 20 km s −1 . It should be noted that some contracting objects are hidden by our choice to present the median in each bin, but even if we were to color-code each cell by its minimum velocity, we would still only be left with a handful of contracting models. We therefore predict that almost all contracting objects are unobservable in the optical lines depicted here due to their high extinction and hence low luminosity.

THE FAILURE OF SIMPLE ABUNDANCE DIAGNOSTICS
The idea of using emission lines from H ii regions to constrain the gas-phase metallicity of the ISM is a very old one (see e.g. Searle 1971, for an influential early paper). The current "gold standard" methods for doing so involve the measurement of the electron temperature Te using line ratios such as the ratio of the [O iii] 4363Å auroral line to the [O iii] 5007Å nebular line that are highly sensitive to the electron temperature. Unfortunately, the auroral lines on which these methods are based are weak and are often not detected in observations of extragalactic H ii regions. Consequently, there is ongoing interest in using the ratio of strong lines, such as [N ii]/Hα, as metallicity tracers (see e.g. Pettini & Pagel 2004). A basic assumption underlying this approach is that variations in the line ratio being used as a metallicity indicator are driven primarily by variations in the relative gas phase abundances of the elements, or in other words that one can cleanly separate the effects of changes in the illuminating radiation field from changes in the elemental abundances. In the case of the [N ii]/Hα ratio, this is often justified by the argument that the similar ionization potentials of N and H make the ratio of N + /H + insensitive to secondary parameters. What these simple assumptions do not account for is the covariance of many secondary parameters in a way that causes significant changes in the [N ii]/Hα ratio at fixed metal abundance. An in-depth examination of the extent to which this covariance confounds commonly-used abundance diagnostics is beyond the scope of our current paper. Nevertheless, the results presented here are already enough to allow us to draw some simple qualitative conclusions.
To begin with, if we look at Figure 9, several points become clear right away. First, there is a good correlation between the age of the youngest cluster in the cloud and the position of the cloud in the BPT diagram. Clouds with young clusters (corresponding to high ionization parameters at the inner shell boundary) have high values of [O iii]/Hβ whereas clouds with older clusters show a strong decrease in [O iii]/Hβ. However, clusters older than ∼ 20 Myr, also populate the region close to the bend in the Kauffmann demarcation line (around log [O iii]/Hβ ≈ 0.2, log [N ii]/Hα ≈ −0.5), but are too faint to be observable and are thus not shown in Fig. 9. We also see that there is a very good correlation between the ionization parameter itself and the lo- Figure 9. Same as Fig. 8 but color coded by instantaneous cloud conditions: (a) the age of the youngest cluster that has formed, (b) the ionization parameter U measured at the inner shell radius, (c) the shell radius, (d) the shell velocity (objects with negative velocities are contracting), (e) the total SFE of all starbursts up to that age, and (f) the derived average star formation rate. cation of the Hii region in the BPT diagram. This result is unsurprising: we know already from static models that there should be a close link between the ionization parameter describing a given Hii region and the line emission from that region. However, our models demonstrate that there is a significant degree of scatter in this relationship, particularly at high [O iii]/Hβ. Observationally, when such scatter is seen, it is often interpreted as a metallicity effect, as changes in the metallicity will generally have the effect of changing [N ii]/Hα while leaving [O iii]/Hβ approximately constant. However, all of our models were calculated at the same metallicity, and so our results clearly demonstrate that metallicity variations are not the only possible cause of the observed scatter.
The top two panels in Figure 9 also illustrate another important result. They demonstrate that there is a good correlation between the value of the ionization parameter at the inner edge of the Hii region and the age of the youngest cluster within the Hii region. Since the cluster age controls the hardness of its spectrum (see Section 2.4), another way of stating this is that the value of ionization parameter is correlated with the hardness of the spectrum. Therefore, we cannot vary these quantities independently when modelling real Hii regions (see also Dopita et al. 2006).

EMBEDDED OBJECTS
As we have shown in Rahner et al. (2017), the minimum star formation efficiency necessary in order for stellar feedback to be able to disrupt a cloud is a function of the cloud mass, metallicity and natal/ambient density. One reason for this is that for a stellar cluster with a fully sampled IMF, the energy and momentum available to drive feedback scale linearly with the mass of the cluster, as M * = SFMcl, but the binding energy of the cloud scales as M 2 cl . Therefore, the amount of feedback per unit of binding energy is not constant, unlike the mass to light ratio of the cluster or the feedback per unit mass ratio, which are both constant for constant SF. Instead, the feedback per unit of binding energy decreases as M −1 cl , leading to larger clouds being more resistant to feedback that smaller clouds.
Furthermore, the natal cloud density also plays an important role. The gravitational binding energy at constant mass is inversely proportional to density, implying that, everything else being equal, denser clouds are harder to destroy. Denser clouds also cool faster, which reduces the effectiveness of stellar wind feedback at early times. On the other hand, they can also be more effective at capturing the full radiation output of the stellar cluster, and so radiation pressure can be more effective in denser clouds. Finally, varying the metallicity also affects the relative importance of winds and radiation pressure. Altogether, the overall impact of feedback depends on the parameters of the natal cloud and the central cluster in a highly complex and non-linear fashion (Rahner et al. 2017;Rahner et al. 2019).
In practice, we find when we vary these parameters that star clusters and their parental clouds form two distinct populations of star-forming regions, with significant differences in both their observational appearance and their star formation histories. We refer to these as unembedded and longterm embedded clusters. Their primary distinction lies in the observability of their stellar clusters. A common observational definition used to distinguish between embedded and unembedded clusters is the mean visual extinction to the cluster: clusters with AV 5.0 are embedded, while those with AV < 5.0 are unembedded (Lada & Lada 2003). We discuss the properties of these two populations in more detail below.

Unembedded clusters
All our natal clusters are initially "embedded" by the AV criterion mentioned above, as expected, but in many cases feedback-driven dynamical evolution rapidly decreases the extinction to the point where we would classify them as unembedded. This is typically associated with the expansion of the shell radius beyond the initial cloud radius, as beyond this point further expansion sweeps up little additional mass from the diffuse ISM but continues to increase the surface area of the shell, resulting in a steady decrease in the column density.
These unembedded clusters have low AV (by definition) and are readily observable because they are still young and hence retain most of the massive stars responsible for ionizing the surrounding gas. The common optical lines used to trace H ii regions (Hα, Hβ, [O iii], [N ii], etc.) are therefore bright. Typically, these models correspond to lower mass star-forming regions with moderate densities (100-200 cm −3 ) that undergo one or two bursts of star formation before successfully disrupting the surrounding cloud and becoming unembedded. Broadly, these regions in which expansion and cloud disruption were successful are the ones explored in most optical surveys of H ii regions.

Embedded clusters
In contrast to the above, we can also identify a class of objects characterized by a long-duration embedded phase, one in which feedback is initially unable to disrupt the cloud and star formation continues over an extended period. This behavior is more likely with more massive clusters, owing to the decrease in the amount of feedback per unit binding energy discussed above. This stands in contrast to the simple picture in which it is often assumed that larger objects lead to monotonically increasing expansion and eventual cloud disruption (see e.g. Murray 2011). Others have also modeled the clouds in massive clusters (Silich & Tenorio-Tagle 2013, including early predictions of failed feedback due to winds (Silich 2004). Advancements were made by including the effect of radiation, but these neglected gravity, or failed to treat radiative coupling in clouds of finite mass with selfconsistently evolving stellar populations (Martínez-González et al. 2014). When considering these effects we find clouds will often not reach Phase III in our evolutionary model (the phase in which the entire cloud is swept up), or if they do, they will undergo re-collapse before expanding far.
We can easily identify these objects from within our set of warpfield models because a significant fraction of massive stars leave the main sequence and die while in the embedded phase (our Phase I and II). Observationally, they are poorly traced by the optical diagnostics that are the main focus of the examples presented in this paper, but they are bright at mid-infrared (MIR) and radio wavelengths. Since warpfield-emp also models the emission at these wavelengths, we can compare the predictions of our models with the results of MIR or longer wavelength surveys, although actually doing so lies outside of the scope of our current paper.
Finally, we note that the existence of embedded clusters that have undergone re-collapse is in itself an important prediction of our model, since it implies that embedded clusters need not be young. Clusters in which feedback has failed to disperse the gas may continue to form stars at a slow rate for an extended period, but may be misidentified as massive young protoclusters.

THE NATURE OF FEEDBACK-REGULATED LINE EMISSION
A natural question that the reader may have at this point is why the approach to modelling H ii region and PDR diagnostics outlined in this paper represents an improvement over the more conventional method of scanning large grids of static models computed by individually varying parameters such as the ionization parameter, the gas density or metallicity, or the age of the central stellar cluster. The central motivation for our dynamical approach is the existence of strong correlations between these different parameters. These correlations mean that in reality we are not free to vary the main parameters controlling the emission independently. For example, the density of the ionized shell is set by the boundary pressure exerted on it by the combined stellar wind from the cluster, which in turn depends on the age, mass and metallicity of the cluster, the current size of the shell, and the phase of the evolution of the region (i.e. whether or not it still retains its hot gas). Therefore, in order for our H ii regions models to be self-consistent, we cannot simply vary the density without also varying some or all of these other parameters. The advantage of warpfield-emp over more static approaches is that it provides us with a simple way of determining which sets of physical parameters are consistent, i.e. which parts of the large parameter space are actually touched upon during the evolution of the feedback-driven shells and which are not. As examples of the power of this approach, we discuss in Section 8.1 the implications of our results for observations of diagnostic lines on the scale of entire galaxies (i.e. observations that do not resolve the individual H ii regions), and in Section 8.2 what warpfield-emp can tell us about the range of cloud-scale star formation efficiencies that we are likely to encounter in reality and what this implies regarding their location in the BPT diagram.
In addition, warpfield-emp also provides us with a means for exploring the impact of particular feedback processes on observable diagnostics of H ii regions since we know at any given time in the models not only how significant the process is at that particular moment, but also how significant it was in the past. This allows us to draw some important conclusions regarding the role played by winds and SNe, which we explore in Sections 8.3 and 8.4 below.

Implications for Mean Galaxy Spectra
Here we note that when weighted by time the distribution of objects tends more toward the mean star-forming galaxy track of Kewley et al. (2013). This effect is particularly apparent when we look at the brightest quartile of the H ii regions (lowest row in Figure 7), which are the objects that generally dominate the mean galaxy spectrum. This is noteworthy, given the wide distribution of the underlying predicted observations. It is a consequence of the self-regulation of star formation in our model clouds. This leads to a covariance between the cluster and cloud parameters which heavily favors observations consistent with this track. A further implication is that this trend is less sensitive to the initial cloud properties, and the emission properties are thus a result of the nature of feedback physics, and the metal abundance of the ISM (which would shift the position of the distribution).

Star Formation Efficiency
In our warpfield models, the star formation efficiency associated with a given burst of star formation is a free parameter. However, because our models account for re-collapse, the total star formation efficiency -i.e. the total mass of stars formed by the end of the life of the cloud -is not completely free. Clouds with very low star formation efficiencies per burst undergo repeated cycles of re-collapse and star formation until they reach a point at which feedback is finally able to disperse them, preventing any further star formation. As a result, our models predict that clouds should have a minimum total star formation efficiency that depends primarily on their surface density, with a value of a few percent for surface densities comparable to local GMCs, increasing to 5-10% in higher surface density clouds . Moreover, although our models do not constrain the maximum star formation efficiency of the clouds, it is reasonable to assume that this will not be much larger than the minimum value required to destroy the cloud. 13 Despite the inherent degeneracy in the limited set of observational diagnostics considered in this paper, our approach leads to some clear predictions about where particular types of cloud should fall in the BPT diagram. In Figure 10, we highlight regions in this diagram populated by unique classes of objects. Objects that have high [O iii]/Hβ but low [N ii]/Hα (blue) correspond to strong starbursts combined with rapidly expanding shells. As Figure 9 makes clear, we expect to find this combination primarily in clouds with high star formation efficiencies. If the simple picture we have sketched above is correct, and clouds rarely develop Figure 10. A generalization of where different classes of objects lie in the observable space. High star formation efficiency, high feedback-to-mass ratio objects (blue) and compact H ii regions and clouds (purple) populate the high [O iii]/Hβ range. Very massive clouds (green) (M cl > 10 6 M ), in which the shell remains highly pressured for an extended period, lie somewhat below the compact regions and tend to have suppressed [N ii]/Hα. Below this, at [O iii]/Hβ −1.0, we find highly extincted H ii regions (red), which also have very low [N ii]/Hα, and extended, old H ii regions that have low ionization parameters (yellow). In dark grey we highlight a region populated by massive, young star-forming regions, with hot interior bubbles. These are deeply embedded super star clusters, and within the realm of observability for deep optical surveys, although distinguishing them from the non-stellar sources that also populate this part of the BPT diagram may prove to be a challenge.
star formation efficiencies much greater than the minimum amount required for cloud destruction, then these objects should be relatively rare in normal spiral galaxies, being confined primarily to regions with high surface density.
The other class of objects that we expect to have high [O iii]/Hβ ratios are small H ii regions exposed to feedback from young clusters (purple). These objects have higher [N ii]/Hα ratios than the strong starbursts and have high [O iii]/Hβ primarily because of their small size and young age. They need not have high star formation efficiencies.
In Figure 8 very massive clouds with M cl > 10 6 M (green) lie somewhat below the small H ii regions in the BPT diagram. They maintain relatively high [O iii]/Hβ values for an extended period, because of the time it takes to completely sweep up the gas in these clouds. Below these lie highly extincted regions, which have low values for both line ratios (red) and extended, old H ii regions with low ionization parameters (yellow), which have low [O iii]/Hβ but somewhat larger [N ii]/Hα. As we have already discussed (see Figures 5 and 7), both of these types of object may be difficult to observe in extra-galactic surveys, particularly if one requires detections of both the [O iii] and [N ii] lines in order to classify them as H ii regions, or when their internal attenuation causes them to become as faint as the diffuse ionized gas (DIG) ).

Significance of Winds
Winds provide an important boundary pressure that plays a major role in determining the density of the feedback-driven shell. This in turn has a substantial influence on how well radiation can couple to the shell. Therefore, even when winds are not a dominant source of feedback in terms of their energy or momentum, they still play a large role in controlling the effectiveness of radiative feedback. Simulations without winds will generally not produce the correct density structure in the gas surrounding the cluster and hence will not correctly capture the effects of radiation pressure.
In terms of observable diagnostics, our models show that the details of the emission produced by the shell are sensitive to the thickness of the shell and hence to the ratio of the wind boundary pressure to the radiation pressure. This changes over time as the cluster ages and the wind-blow bubble cools and expands.
Even before we explore the impact of changing metal abundances, and the change in wind strength that follows from this, it is clear that the change in the relative thickness of the H ii region due to the evolution of the cluster and the expansion of the bubble is a detail that cannot be ignored, as it can result in orders of magnitude changes in the relative strengths of low and high ionization lines.

The importance, or lack thereof, of SNe
An unexpected result of our models is our finding that many of the evolving clouds, particularly the least massive ones, meet our dissolution criteria immediately after SNe start exploding, but that the SNe themselves are not responsible for destroying the cloud. Recall from Section 2.1 that our primary criterion for cloud dissolution is that the densest portion of the cloud must have a density n < 1 cm −3 for at least 1 Myr. In the clouds in question, the feedback-driven shell sweeps up the entire cloud before any of the massive stars in the central cluster reach the end of their lives (i.e. the clouds reach Phase III in the terminology of Section 2). While all of the stars in the cluster remain alive, the density of the gas in the shell remains above 1 cm −3 owing to the confining pressure of the winds and radiation. However, once some of these stars start to explode as SNe, the momentum available in the winds and the radiation that couples to the shell rapidly decreases, despite the positive contribution made by the SNe themselves. As a result, the interior pressure of the cloud drops quickly, resulting in the peak density falling below our dissolution criterion.
Clouds that behave in this fashion will tend to leave behind relatively bright clusters which continue to ionize their surroundings. Moreover, these clusters will be associated with evidence of recent SN explosions. Observationally, one might therefore incorrectly conclude due to the timing that it is the SNe that are responsible for the destruction of the cloud, when in reality, it is the drop in internal pressure due to the aging of the cluster that leads to the shell density dropping and the shell becoming difficult to observe (or indistinguishable from diffuse ionized gas).

SUMMARY
In this paper, we have presented the warpfield-emp framework, which is comprised of warpfield, which we use to model the dynamics of feedback-affected gas around star clusters, the Cloudy spectral synthesis code, which we use to model the line and continuum emission produced by this gas and its associated dust, and polaris, which we use to account for the effects of dust extinction and line opacity on the emergent radiation.
warpfield-emp allows us to model the time evolution of feedback in molecular clouds without ignoring any of the physical processes that are important for regulating emission from the clouds. The list of physical processes included in the model includes stellar winds, radiation, supernovae, gravity, thermal conduction, cooling of the hot wind-blown bubbles directly via radiation, and indirectly by the introduction of cold entrained shell material. All of these processes are interconnected and highly non-linear, and hence cannot be studied reliably in piecemeal fashion.
The main approximation made in the warpfield-emp approach is the assumption of 1D spherical symmetry (although the model includes non-symmetric effects such as shell fragmentation in a subgrid fashion), a common assumption in nearly all extra-galactic photoionization modeling. Making this approximation reduces the computational cost of the model to the point where it becomes plausible to study large samples of clouds, something which will remain out of reach in 3D simulations for the forseeable future. warpfield-emp allows us to study large parts of the cloud and cluster parameter space. And so it is no longer easily justifiable to treat aspects of H ii region physics such as the ionizing SED, ionized gas pressure and density, or the size of the region as independent parameters, nor to ignore the extensive properties of H ii regions and PDRs, such as observed luminosity and cloud mass. As we have shown in this paper, these parameters are not independent and ignoring the substantial covariance between them can easily lead to erroneous conclusions.
We have explored a large parameter space spanning low to high cloud and cluster masses and a range of different cloud densities as a first step toward understanding how nature draws from the cloud and cluster mass function and convolves that into observable diagnostics. We have compared the predicted emission from our model clouds to BPT diagrams of H ii regions in NGC 628 as a first step toward verifying the approach. We find very good agreement between our models and the observational data once we properly account for the observational selection effects. Notably, we reproduce not only the shape of the locus of observed H ii regions in the BPT diagram but also the scatter. This scatter has often been attributed to metallicity variations, but our constant metallicity models demonstrate that it may simply reflect differences between the evolutionary state of different H ii regions.
In a pair of companion papers (including Pellegrini et al. 2019) that use our warpfield-emp framework as the basis for a sophisticated population synthesis treatment of galactic emission, we show that the models also do a good job of reproducing the intensity and structure observed in all-sky maps of Galactic Hα emission as well as measurements of Faraday rotation.
Our models allow us to predict the emission from a large population of evolving clouds in a fully self-consistent fashion. Despite the inherent degeneracy in the diagnostics considered in this paper, they allow us to draw a number of important conclusions: (i) There is a causal relationship between the feedback, and cluster/cloud evolution which determines the ionizing SED, ISM pressure, ionization parameter and cloud structure. The coupling of these parameters is set by a combination of microphysics (which can be modelled from first principles) and the cumulative effect of feedback on the global dynamics of the clouds. Together, these determine both the time evolution of the cloud properties and the star formation history of the clouds.
(ii) Constraining photoionization model parameter spaces with dynamics breaks key degeneracies in very simple models between age, size and luminosity.
Replacing the old way of modeling nebulae (e.g. by computing static grids of quantities such as ionization parameter, metal abundance, stellar effective temperature, etc.; see e.g. Kewley et al. 2001;Kewley & Dopita 2002;Kewley et al. 2013;Pellegrini et al. 2011;Byler et al. 2017) with actual clouds that evolve self-consistently breaks some important degeneracies (although it also introduces new ones related to natal cloud properties). For expected cloud properties, even at constant metal abundance, photoionization tracks of realistic objects are not clean, as observed in the common [N ii]/Hα vs. [O iii]/Hβ plane.
(iii) In massive clouds and high pressure regions, collapsing objects are often unobservable at optical wavelengths owing to their high extinctions. Surveys which identify H ii regions using common optical lines (e.g. Hα) will tend to miss regions undergoing re-collapse, an observational bias which is important to account for when using H ii region properties to derive feedback properties. A corollary of this is that embedded clusters observed at mid-or far-infrared wavelengths need not be young, but could in some cases contain stellar populations with ages up to 30 Myr.
(iv) Deeply embedded super star clusters appear as faint AGN.
The location of deeply embedded massive clusters in the BPT diagnostic (see Figure 10) presents challenging issues to observers who wish for a simple diagnostic to separate active galactic nuclei (AGN) from star formation activity. As shown, the large amount of differential extinction, combined with large internal gradients in ionization parameter can result in star-forming clouds having diagnostics consistent with faint AGN, leading to an unknown number of missing/or wrongly classified super star clusters (SSCs).
(v) Simple line diagnostics, like [N II]/Hα, have substantial scatter at constant ionization parameter and abundance.
The primary cause of this scatter is changes in the boundary pressure between different regions (or within the same H ii region at different times), which have a substantial effect on the thickness and density of the ionized shell. Pressure evolution causes H ii regions in clouds of finite mass to migrate between being density-bounded and being radiation-bounded (i.e. between having a low or a high optical depth to ionizing radiation).
(vi) Clouds with low star formation efficiencies may undergo multiple bursts of star formation.
If the star formation efficiency is too low, feedback from the star cluster may be unable to unbind the cloud, which will eventually re-collapse and form new stars (see also Rahner et al. 2019). This implies that some "rejuvenated" clusters will contain multiple stellar populations with ages separated by 1-10 Myr. A striking example of this may be 30 Doradus in the Large Magellanic Cloud (see e.g. Rahner et al. 2018). Below 1 µm, the emission from these rejuvenated clusters is dominated by their youngest population, hiding older generations of stars in optical and UV observations. This has serious but largely unexplored implications when measuring star formation rates and efficiencies from clusters of multiple ages, with mass traced by different methods (e.g. Hα vs Vor B-band magnitudes).
when ionizing radiation pressure would otherwise drive gas densities high in the cold gas layer. We treat turbulent pressure support using a similar approximation. We assume that it does not significantly alter the dynamical evolution of the shell, but allow the user to account for small-scale turbulent pressure support in the shell. We do not include this effect by default because although molecular clouds are ubiquitously observed to be turbulent, the degree to which this turbulence acts to provide smallscale support is highly unclear (see e.g. the discussion in Mac Low & Klessen 2004).
Our RT and line emission calculations begin at the boundary of the wind blown bubble and the H ii region shell. In the absence of magnetic fields and turbulent support, we assume that the thermal pressure here is equal to the boundary pressure of the bubble P b : However, when turbulence and a magnetic field are also included, we instead equate the boundary pressure with the total shell pressure: The assumption of energy equipartition between the magnetic field and the turbulent field implies that and if we use a scaling law for the magnetic field in the form of with B0 = 48.7 µG, n0 = 116.1 cm −3 , γmag = 4/3, and Ti = 10 4 K. In calculations with B-fields and turbulence we use this equation, which can be numerically solved for the ion number density at the inner edge of the shell n(R), instead of Eq. B1.

APPENDIX C: SUPPLEMENTAL MATERIAL
The following supplemental material provides some additional examples of evolutionary tracks along the lines of those shown in Figure 5, but for different cloud and cluster properties. In Figure C1 we show an example of a low star formation efficiency cloud, with a high density. This results in many recollapses, as the cloud is not destroyed by feedback after many Myr. In Figure C2 we show a cloud with moderate density of 300 cm −3 , with a high 10% SFE.
Here feedback is sufficient to drive to cloud to large radii of almost 100 pc. The cloud becomes optically thin to ionizing radiation and the [N ii]/Hα ratio decreases rapidly after shell expands beyond the natal cloud boundary. Finally, SNe mark a rapid decrease in pressure. While they do not destroy the cloud, which has already been disrupted by winds and radiation, the pressure loss associated with them leads to the density of the shell falling below that of the warm ionized medium. At this point, we consider the shell to have been destroyed.

APPENDIX D: DEFAULT EMISSION PREDICTIONS
The following is a list of the default set of line and continuum predictions made in warpfield-emp. These include line and continuum emission, as well as flux-weighted opacities from dust and atomic and molecular processes. These can be calculated integrated across the nebula, or in projection. Note that this set does not include all of the lines that warpfield-emp is capable of modelling; instead, it represents the subset that we expect to be most useful as tracers of H ii region, PDR or cloud properties. Figure C1. Same as Figure 5 but for a natal density of 500 cm −3 , cloud mass log 10 M cl = 5.50 and SF =1%.