{'fraction_non_alphanumeric': 0.023086984957488555, 'fraction_numerical': 0.017069980379332897, 'mean_word_length': 4.759794073329985, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the relative importance of different microphysics on the D-type expansion of galactic HII regions

Radiation hydrodynamics (RHD) simulations are used to study many astrophysical phenomena, however they require the use of simplified radiation transport and thermal prescriptions to reduce computational cost. In this paper we present a systematic study of the importance of microphysical processes in RHD simulations using the example of D-type HII region expansion. We compare the simplest hydrogen-only models with those that include: ionisation of H, He, C, N, O, S and Ne, different gas metallicity, non-LTE metal line blanketed stellar spectral models of varying metallicity, radiation pressure, dust and treatment of photodissociation regions. Each of these processes are explicitly treated using modern numerical methods rather than parameterisation. In line with expectations, changes due to microphysics in either the effective number of ionising photons or the thermal structure of the gas lead to differences in D-type expansion. In general we find that more realistic calculations lead to the onset of D-type expansion at smaller radii and a slower subsequent expansion. Simulations of star forming regions using simplified microphysics are therefore likely overestimating the strength of radiative feedback. We find that both variations in gas metallicity and the inclusion of dust can affect the ionisation front evolution at the 10-20 per cent level over 500kyr, which could substantially modify the results of simplified 3D models including feedback. Stellar metallicity, radiation pressure and the inclusion of photodissociation regions are all less significant effects at the 1 per cent level or less, rendering them of minor importance in the modelling the dynamical evolution of HII regions.


INTRODUCTION
Radiative feedback has the potential to influence the morphological evolution of star forming regions and to induce or inhibit star formation (for a recent review see Dale 2015). The radiation field from massive stars ionizes the surrounding gas, heating it and causing expansion of the hot high pressure region (Spitzer 1978;Hosokawa & Inutsuka 2006;Raga et al. 2012). This can result in the dispersal of material, hindering the formation of stars, or can collect/destabilize material and potentially trigger star formation. A large number of factors contribute to the effectiveness of radiative feedback (making models difficult, because there is a huge E-mail: thaworth@ast.cam.ac.uk parameter space of initial conditions and microphysical complexity) and feedback processes are expected to take place over timescales of order 1 Myr in systems of complex geometry (making observations difficult, due to projection effects and the challenge of constructing a dynamical picture from a single snapshot). Consequently the impact of radiative feedback on star formation is still not clear, particularly in any quantitative sense (see e.g. Dale et al. 2015).
Radiation hydrodynamics (RHD) calculations have been used to investigate the effect of radiative feedback in a host of "star formation" scenarios. The external irradiation of an isolated cloud (e.g. Esquivel  and collapse (e.g. Dale et al. 2007) and the radiatively driven evolution of turbulent media (e.g. Arthur et al. 2011;Walch et al. 2012; have been studied using numerical models, providing a phenomenological picture of the impact of radiative feedback on the evolution of both the gaseous and stellar content of star forming regions. These RHD models all necessarily use (to varying extents) a simplified treatment of radiation transport/photoionisation due to the complexity of the microphysics and finite computational resources available. In recent years some effort has been invested into understanding how important these assumptions are. For example  investigated the effect of including polychromatic and diffuse field radiation in models of the radiatively driven implosion of clouds, finding that the diffuse radiation field can significantly modify the results.  also found that the assumption of photoionisation equilibrium can affect the results of RHD calculations of the external irradiation of a turbulent medium. They found that non-equilibrium photoionisation was required to detatch the tips of elephant trunks to form Bok globules (Bok & Reilly 1947). More recently Sales et al. (2014) investigated the effect of radiation pressure, finding that it is a secondary effect compared to photoionisation. Geen et al. (2015) also studied the relative impacts of winds, ionising radiation and supernova feedback in a series of 1D models, which they used to summarise the energetic feedback into the ISM from a 15 M star.
There are further approximations that have not been formally investigated, for example the assumption that the gas is hydrogen-only (thus neglecting cooling from forbidden line transitions and heating and cooling from helium) and using a simplified thermal balance that calculates the temperature as a function of hydrogen ionisation fraction (though e.g. Raga et al. 1999;Mellema et al. 2006;Mackey & Lim 2010;Miao et al. 2006;Tremblin et al. 2012, com-pare the heating and cooling rates). There are also variations in the gas or stellar metallicity, dust and photodissociation regions/FUV heating. The impact of these approximations must be investigated to understand by how much and why simple models differ from more complex ones.
We can easily illustrate how different approximations might be expected to affect the evolution of H ii regions by considering the classic system of a massive star at the centre of a uniform density medium. The star rapidly ionises a sphere of the surrounding gas. If the gas is hydrogen only and we invoke the on-the-spot (OTS) approximation, under which the diffuse radiation field is neglected, then the radius of this initial bubble is the Strömgren radius where N ly , ne and αB are the number of ionising photons emitted by the star, gas electron density and the case B recombination coefficient (for recombinations into all states other than the ground). The subsequent D-type expansion of the H ii region was described by Spitzer (1978) as where cI is the sound speed in the ionised gas. If we consider a case B recombination coefficient sensitive to temperature as T −0.8 (see equation 5) then the Strömgren radius varies as T 0.27 . This will affect the expansion rate at times when 7 c I t 4 rs 1. At times when 7 c I t 4 rs >> 1 we expect that the expansion of an H ii region will vary with the temperature in the ionised gas as rI ∝ T 2/7 I t 4/7 and the number of ionising photons as rI ∝ N 1/7 ly t 4/7 . Although these dependencies are quite weak, departures from simple estimates in these quantities may add up over time to give substantial H ii region expansion differences. Furthermore, if we improve our model by considering gas that is not hydrogen only and includes the diffuse field then both the recombination coefficient and the electron density will change, further modifying the Strömgren radius which affects the expansion. There are hence multiple factors that could affect the expansion of H ii regions when improving the microphysical treatment.
In this paper we aim to explore how different stellar and gas metallicities, radiation pressure, dust and photodissocation regions affect the D-type expansion of H ii regions, which we interpret in terms of the simple analytic expectations mentioned above. Specifically we will investigate the effects of these processes in galactic H ii regions rather than those studied in the epoch of reionisation and make the distinction based on the gas densities, gas constituents and the ionising fluxes from the sources considered.

NUMERICAL METHOD
We use the Monte Carlo radiation transport and hydrodynamics code torus (Harries 2000;Kurosawa et al. 2004;Rundle et al. 2010; to perform the calculations in this paper. The torus RHD algorithm uses operator splitting to separate grid-based hydrodynamics and Monte Carlo photoionisation ). The primary advantage of this method is that all of the features available to a dedicated Monte Carlo radiation transport code are available in RHD calculations. The disadvantage is that this approach is computationally expensive, however it can be efficiently parallelized using a range of techniques (Harries 2015). Details and testing of the RHD algorithm are provided in  however we include a summary here for completeness since this paper predominantly explores different physical processes that the code can include. We also use the coupled torus-3dpdr code, which is discussed in detail in Bisbas et al. (2015) and also summarised here.

Hydrodynamics
We use a flux conserving, finite difference hydrodynamics algorithm. It is total variation diminishing (TVD) and makes use of the van Leer flux limiter (van Leer 1979) and a Rhie-Chow interpolation scheme to prevent odd-even decoupling (Rhie & Chow 1983). torus is capable of treating point source (Harries 2015) and self-gravity, the latter of which is calculated using a multigrid method, though we do not include gravity in the models in this paper.

Photoionisation and thermal balance
torus uses a photoionisation scheme similar to that of Ercolano et al. (2003) and Wood et al. (2004) which in turn are based on the method presented by Lucy (1999). Packets of photons at constant frequency and that carry constant energy (but whose members vary in number with frequency) are propagated throughout the computational grid. As they traverse grid cells they trace a path length l and modify the time-averaged radiation energy density U in the cell by where V is the cell volume, c is the speed of light, Jν is the mean intensity and ∆t is the time over which the averaging takes place. The update to the radiation energy density takes place following any photon packet 'event' which, as well as absorption, includes traversal of a cell boundary. Following an absorption the photon packet is re-emitted with random frequency and direction under the principle of detailed balance, continuing a random walk through the grid until it escapes. The spectrum for diffuse field photons is constructed using 1000 frequencies that strategically sample Lyman continuum and helium ionising photons, hydrogen recombination lines, forbidden lines and the dust continuum (the diffuse field is thus dependent on the species included and ionisation and temperature state of the gas). Once all photon packets in the calculation have escaped the radiation energy density is used to solve the ionisation balance equation (Osterbrock 1989) which, in terms of Monte Carlo estimators, is given by where n(X i ), α(X i ), aν (X i ) and ne are the electron number density, recombination coefficient and absorption cross section of ion X i and the electron density respectively. The approach of using the crossing of cell boundaries as a photon packet event has the advantage that photon packets contribute to the estimate of the radiation field without having to undergo absorption events, thus even very optically thin regions are properly sampled. For the simplest models in this paper we assume that the gas is either entirely atomic or ionised hydrogen. The radiation field is monochromatic and we use the OTS approximation, including the same case B recombination coefficient as used by Bisbas et al. (2009) Our simple models employ a common (e.g. Bisbas et al. 2011) simplified thermal balance calculation where the temperature in cell j is a twotemperature interpolation function of the ionisation fraction of atomic hydrogen ηj Tj = Tn + ηj(Tio − Tn) where Tn is the prescribed temperature of fully neutral gas and Tio the prescribed temperature of fully ionized gas (10 and 10 4 K respectively). We retain the assumption of photoionisation equilibrium in all models. In contrast to the simple calculations, in the detailed photoionisation models we include a range of atomic constituents: hydrogen, helium, carbon, nitrogen, oxygen, neon and sulphur, for which we solve the ionisation balance using equation 4. The levels that we treat for metals are C (I-IV), N (I-IV), O (I-III), Ne (II-III), S (II-IV). The hydrogen, helium and C IV recombination rates used by torus are calculated based on . Other radiative recombination rates are calculated using fits to the results of Nussbaumer & Storey (1983), Pequignot et al. (1991) or Shull & van Steenberg (1982). The photoionisation cross sections of all atomic species in this paper are calculated using the phfit2 routine from . We note that we assume rather than explicitly calculate the abundance of helium and metals and use this assumed abundance to calculate the ionisation structure. The detailed photoionisation model temperature is calculated by finding the temperature at which the heating and cooling rates in each cell match. The gas heating rate from hydrogen and helium in a given cell is calculated based on the sum of trajectories of photon packets through the cell (Wood et al. 2004). The dust heating is calculated separately, but using the same method (Lucy 1999).
The cooling rate is initially calculated for the maximum and minimum allowed temperatures in the calculation (3×10 4 K and 10 K respectively by default in torus). This is then refined by bisection until the cooling rate matches the heating rate. For gas, the cooling processes considered are those from free-free radiation, hydrogen and helium recombination and collisional excitation of hydrogen and metals. For dust, there is blackbody radiative cooling. We iterate over the ionisation and thermal balance calculations until the fractional change in both the ionisation fraction and temperature is less than 10 −2 . Gas and dust are thermally (but not dynamically) decoupled, having their thermal balance solved independently.
In the detailed calculations that are not investigating the effect of dust we set the dust to gas ratio to the negligibly low value of 10 −20 .

OSTAR2002 Spectral models from TLUSTY
For stellar spectral models, we use the non-LTE, metal line blanketed, plane parallel radiation transport and hydrostatic equilibrium "OSTAR2002" models of Lanz & Hubeny (2003), calculated using the code tlusty. We use three sets of grids for stellar metallicities of Z = 0.5, 1 and 2 Z . Each grid consists of 69 models, spanning temperatures from 27500 K to 55000 K and surface gravities from 3.0 log(g) 4.75. The actual spectrum used in a calculation is derived by interpolation between the two grid spectra with properties closest to that of the star in our model. In Figure 1 we show an example blackbody spectum and that interpolated from OSTAR2002 for a star at 45000 K, a radius of 10.9 R and various stellar metallicities. The upper panel shows the spectrum over a large wavelength range and the lower panel shows the spectrum from the Lyman limit.
2.3.1 Frequency sampling of the stellar and diffuse field spectra torus converts the source spectrum into a cumulative probability distribution as a function of frequency. Generating a random number in the range 0:1 and mapping this onto the cumulative probability distribution then yields a frequency with probability appropriate to the source spectrum. The tlusty spectra consist of 193434 frequencies per model and the blackbody spectra consist of 1000 frequencies. The blackbody spectrum is logarithmically sampled from 10 to 10 7Å .

Radiation pressure
Harries (2015) discusses and tests the treatment of radiation pressure by torus in detail. In summary the radiation pressure force is calculated using Monte Carlo estimators. This estimate is calculated as photon packets are propagated over the grid in the photoionisation component of the calculation and so if already doing the photoionisation, is essentially obtained for no additional computational cost. This technique accounts for polychromatic radiation and anisotropic scattering and works in both the free-streaming and optically thick regimes. The calculated radiation pressure appears as an additional force term in the hydrodynamic component of the calculation.

Dust
Historically, torus has been used to study discs around young stars where the gas and dust can be assumed to be thermally coupled. For modelling of H ii regions, which are at much lower densities, we must thermally decouple the dust from the gas. This has been done by other codes (e.g. Ercolano et al. 2005;Pavlyuchenkov et al. 2013), but by torus for the first time in this paper. We assume spherical silicate dust particles that follow a standard interstellar medium power law size distribution (e.g. Mathis et al. 1977) of the form where the number of grains of radius a is where c and q are constants. The dust optical constants are taken from Draine & Lee (1984). We use a pre-tabulated Mie-scattering phase matrix. At present our dust treatment does not include photoelectric heating or resonant line transfer. We assess the impact of this approximation by comparing with the more advanced (in terms of photoionisation) Monte Carlo photoionization code mocassin (Ercolano et al. 2003(Ercolano et al. , 2005 which does include these dust processes.
The model that we use for testing is the 1D HII40 Lexington benchmark (see Ferland 1995;Ercolano et al. 2003;, only with the inclusion of dust. The dust to gas mass ratio is 10 −2 and we use Draine silicate grains. The density is 100 mH cm −3 and the metal abundances used by both codes are as given in Table 1. A comparison of the gas temperature distribution (all that really matters for these dynamic calculations, as opposed to the line intensities) as computed by torus and moccasin is given in Figure 2. Clearly the thermally decoupled gas and dust model from torus calculates a very similar H ii region radius and temperature to mocassin. The extent of the ionised gas has been reduced compared to the simulation in which there is no dust. Beyond the ionisation front there is slight heating (of order tens of Kelvin) in the neutral gas by torus, but this is weaker than in the mocassin calculation, where there neutral gas is heated to ∼ 300 K. We do not expect this downstream heating to have much effect on the dynamics compared to the more dramatic effect on the ionisation front radius. We note that this additional heating in the mocassin calculation is very similar to photodissociation region heating, which will be investigated in this paper, so we will be able to gauge its impact.

TORUS-3DPDR
We recently coupled torus with the 3d-pdr code of Bisbas et al. (2012), the details and testing of which are given thoroughly in Bisbas et al. (2015). This approach is novel since the UV field used in the PDR calculation is determined explicitly by the Monte Carlo radiation transport (typically the Draine field is a free parameter for PDR codes). The coupled code is fully integrated, with the PDR calculations taking place on the torus grid. This means that we can calculate the ionisation state of multiple species in the H ii region as well as the chemical and thermal structure of the PDR and then use the resulting thermal properties as pres-   Motoyama et al. (2015) recently presented a 2D hydrochemical code capable of modelling PDR chemistry and evolving the gas dynamics. The difference with our code is that we can directly calculate the UV field over all space and can also calculate the properties in the H ii region. Unfortunately such a calculation is computationally very expensive. torus-3dpdr uses the most recent UMIST 2012 chemical network database, consisting of 215 species and more than 3000 reactions (McElroy et al. 2013). However at present we use a reduced network of 33 species and 330 reactions, primarily to reduce computational cost. Nevertheless, RHD simulations including even this reduced network of species and reactions, in conjunction with photoionisation, are novel. Note that we assume equilibrium in both the H ii region and the PDR.
Modelling of the PDR should give rise to higher temper-atures (a few hundred Kelvin rather than 10) at the boundary of the H ii region. This will allow us to investigate the effect of heating seen by dust just beyond the ionisation front by mocassin that is not replicated by torus (see Figure 2).

TEST MODEL SPECIFICATIONS
The model that we use as a testbed is the classic expansion of an H ii region about an ionising star in a uniform density medium (Spitzer 1978;Dyson & Williams 1980;Hosokawa & Inutsuka 2006;Raga et al. 2012;Bisbas et al. 2015). This system is 1D spherically symmetric, reducing computational cost and thus allowing us to incorporate a large range of microphysics in this paper. Specifically we only consider the early phase expansion (Bisbas et al. 2015).
We consider a star of effective temperature 45000 K, radius 10.9 R and mass 63.8 M . These stellar parameters are taken from Diaz-Miller et al. (1998) and is in the regime where non-LTE model atmospheres are required. For a model star with these parameters we consider the expansion of the H ii region into two different ambient densities of 100 and 500 mH cm −3 which we refer to as low and high density respectively. For this scenario we run a host of models with different physical prescriptions: simple photoionisation and full photoionisation (see section 2.2), non-LTE spectral models, different gas metallicities, radiation pressure, dust and treatment of photodissociation regions. We study the gas and stellar metallicity effects separately so that we can determine which, if either, the dynamics is most sensitive to. In reality the gas and stellar metallicities will not be independent. A summary of the physical parameters in the model is given in Table 1.

RESULTS AND DISCUSSION
We compare the models by tracking the ionisation front position (defined as the point at which the hydrogen ionisation fraction is a half) as a function of time, this is shown for all low and high density models in Figures 3 and 4. We now discuss each set of approximations/physical processes in turn.

Simple versus full photoionisation
The simple photoionisation scheme assumes an ionised gas temperature of 10 4 K. If the temperature resulting from the detailed photoionisation differs from this then the expansion rate of the H ii region will differ, since the Strömgren radius varies as rs ∝ T 0.27 (for a case B recombination coefficient) and at later times during the expansion varies as rI ∝ T 2/7 (see equation 2). The simple scheme also considers hydrogen only gas, which means that the electron density compared to a calculation with helium and metals may also differ. Finally the simplified scheme also uses the OTS approximation and so there may be differences between the case-B recombination coefficient used (see equation 5) and the net effective recombination coefficient for the calculation with multiple species.
The top left panels of Figures 3 and 4 show the ionisation front position as a function of time for the Spitzer analytic solution (equation 2) as well as for our simplified and  full photoionisation calculations with torus for the 100 and 500 mH cm −3 density calculations. We also include the analytic expression describing the expansion given by Hosokawa & Inutsuka (2006) (labelled H-I), which is derived by solving the equation of motion of the shell of material swept up by the expanding H ii region which, under the assumption that Po Pi, results in The difference between the H-I and Spitzer solutions is that the former solves the equation of motion of the shell (thin shell approximation) and the latter considers pressure balance between the H ii region and ambient medium at each point in time).

Initial Strömgren radii
Looking at the top left panels of Figures 3 and 4, the simplified models agree with the Strömgren solution for the position of the onset of D-type expansion in both density regimes, however for the detailed model this radius is smaller by 7.2 and 6.5 per cent in the 100 and 500 mH cm −3 models respectively. Given that the number of ionising photons is the same in each model, the difference can only come from a difference in the electron density or recombination rate from equation 1. We compared the electron density for the two models at the onset of D-type expansion and found they are very consistent. The difference must therefore be arising from a higher effective recombination rate across all species relative to the case B recombination coefficient used to calculate the Strömgren radius analytically. The recombination rate (equation 5) is temperature dependent, so the difference could be explained if the simple and detailed photoionisation schemes result in different gas temperatures. We plot the gas temperature as a function of radius from the ionising source in Figure 5. The simplified calculation results in higher average gas temperatures, implying a lower recombination rate and hence a larger Strömgren radius (as we observe in our models).
In Figure 6, which we will discuss more below, we compare our simplified and detailed photoionisation results against analytic solutions using gas temperatures representative of the model average in the ionised gas. When we do this the initial I-front position of the detailed photoionisation models is in good agreement with the Strömgren equation. So we can conclude that the difference in Strömgren radius is explained by a difference in the gas temperature and hence recombination rate.
We note that the Hii region radius calculated by torus for detailed photoionisation models has been shown to agree with other photoionisation codes, e.g.  where we compare with the cloudy code (Ferland et al. 1998).

Expansion rates
We remind the reader that the top left panels of Figures 3  and 4 show the ionisation front position as a function of time of the models being discussed presently. The StarBench 1 code comparison project D-type expansion test (which uses the simplified photoionization physics and towards which torus contributed results) found that numerical codes agree with the Spitzer solution early on, but eventually overtake it to agree more closely with the Hosokawa-Inutsuka solution. This behaviour is reproduced in the simple photoionization models in this paper, though the expansion is slightly slower than expected in the lower density regime. We expect that this is due to the comparatively low resolution (256 cells) used for the simulations in this paper compared to the 1024 cells used in the 1D simulations of Bisbas et al. (2015). Although we could increase the resolution for the simplified calculations (indeed our contributions to the StarBench Dtype test were at higher resolution) this would make consistent resolution calculations including the PDR very computationally expensive. The effect of moving to full photoionisation is not only to reduce the radius of D-type onset (see section 4.1.1), but also to reduce the expansion rate of the H ii region. In Figure 5 we show the temperature of the detailed and simple models just before the onset of D-type expansion and after 250 kyr. Throughout most of the H ii region the temperature in the ionised gas of the full photoionisation model is around 2000 K lower than that of the simplified model, except for the temperature peak near to the ionisation front. These cooler average temperatures are due to the efficient forbidden line cooling, which is only included in the detailed model. The peak in temperature close to the ionisation front arises where coolants with higher ionisation potentials recombine. The simple and detailed models have average gas temperature of 10 4 and 8000 K respectively. In Figure 6 we use these average gas temperatures in the Spitzer solution and compare with our numerical results. This demonstrates that differences in the D-type expansion between the simple and detailed photoionisation models are simply explained by the differences in the ionised gas temperature.

Stellar spectral models
Although the different model spectra shown in Figure 1 appear significantly different, the integrated ionising flux over the spectrum is actually very similar in each case for stars of this effective temperature (see table 2 in this paper and Figure 15 from Lanz & Hubeny 2003). The expansion of H ii regions about stars with these different spectra is therefore very similar, as we show in the top right panels of Figures 3  and 4. The difference in the initial Strömgren radius is negligible and after 500 kyr of expansion the difference in H ii region radius is only of order 1-3 per cent in both density regimes.
The lower panel of Figure 1 shows that harder photons are produced in the tlusty models, however this has very little effect on the ionisation state or dynamical evolution, only slightly broadening the ionisation front. This is in agreement with what was found by  in 3D models of radiatively driven implosion, where the effect of including harder radiation was isolated. Although hard radiation can have important consequences at galactic or cosmological scales (e.g. Shapiro et al. 2004;Iliev et al. 2006), on scales of up to tens of parsecs we find that its effect is negligible.

Gas metallicity
The gas metallicity determines the amount of metal line cooling and therefore the temperature in the ionised gas. Since the Strömgren radius varies as T 0.27 and at later times rI ∝ T 2/7 the expansion rate of H ii regions is expected to Figure 6. Demonstrating that the differences between the simplified and detailed photoionisation models can be explained by their different ionised gas temperatures. The upper and lower panels show results for the 100 and 500 m H cm −3 density models respectively. The analytic solutions are just the Spitzer equation with an ionised gas temperature of 10000 K and 8000 K for the simple and detailed photoionisation models respectively. be metallicity dependent. We find that this is the case in our models, which we show in the middle left panels of Figures 3 and 4. Included are calculations at the base metallicity (which we call zG = 1, that is the same as that used in the photoionisation HII40 Lexington benchmark) as well as zG = 0.5 (representative of earlier universe, LMC) and zG = 2 (representative of the Galactic centre). The lower metallicity expansion is much faster than the higher metallicity, giving a zG = 0.5, 500 kyr H ii region radius larger than that for zG = 2 by 18 and 16 per cent in the 100 and 500 mH cm −3 models respectively.
In Figure 7 we show the radial temperature profile of the different metallicity H ii regions at a time close to the onset of D-type expansion. The temperature in the ionised gas and the Strömgren radius are both decreasing as a function of gas metallicity. This is because at lower metallicity there is weaker metal line cooling. Figure 8 shows the ionisation front position as a function of time for the models of different gas metallicity, as well as analytic plots which are generated from the Spitzer solution and Strömgren radius, only setting the ionised gas temperature to the average seen in Figure 7. The agreement is good enough for us to conclude that the differences in H ii region expansion rate with metallicity are dominated by differences in the ionised gas temperature.
Since there is good agreement between the theoretical expressions and simulations when the correct temperature is used, we fit our temperature as a function of metallicity to obtain an approximate, simplified thermal calculation that accounts for the gas metallicity, improving on the simple prescription given by equation 6 T = Tn + 1.1 × 10 4 − 3.8 × 10 3 (z/zo − 0.5) 0.839 − Tn η (10) where z/zo is the gas metallicity relative to that of the Lexington benchmark (a standard Milky Way star forming region abundance), Tn is the neutral gas temperature and η is the hydrogen ionisation fraction. An example potential application of this approximate simplified temperature calculation would be to use it 3D models such as those of Dale et al. (2013) or Walch et al. (2013) (both of which use equation 6) to study how radiative feedback affects molecular clouds in regions of different metallicity, or as a function of cosmic time.
In reality the gas and stellar metallicities will not be independent. We investigate them separately in this paper to identify which of the two components dominates any differences in the H ii region evolution. Based on our simulations we find that feedback dynamics are much more sensitive to the gas metallicty than changes in the stellar metallicity.

Radiation pressure
The ionisation front evolution for models that include radiation pressure is given in the middle right panels of Figures  3 and 4. The radiation pressure simulations are just the "detailed" models with the addition of the radiation pressure formulation developed by Harries (2015). In both density regimes the effect of including radiation pressure is essentially zero (<1 per cent). This is in agreement with Sales et al. (2014), who find that radiation pressure effects are secondary to photoionisation in D-type expansion models. Further note that Raga (2015) finds that hundreds of O stars are required for radiation pressure to become important to the system dynamics, such as in the Tarantula nebula in the LMC.

Dust
Dust absorption reduces the ionising photon budget and hence the size of the Strömgren radius (see equation 1). In Figure 2 we also saw that it increases the temperature downstream of the ionisation front in the mocassin results (though this is missing from our treatment of the dust in this paper, the consequences of this downstream heating effect will be gauged in our discussion of PDR heating in section 4.6). The ionisation front expansion for the dust models is given in the bottom left panels of Figures 3 and 4. As we found when testing in Figure 2, inclusion of dust leads to the initial Strömgren radius being reduced. This affects the expansion rate of the H ii region at times when the Strömgren term is important, i.e. when 7 c I t 4 rs 1. Although the ionising photon budget is reduced, the gas temperature remains very similar (see Figure 2) and so at later times when rI ∝ T 2/7 the expansion behaves in a similar manner either with or without dust. The effect of dust, minus any effect from downstream heating, is therefore simply to reduce the radius at which D-type expansion begins, resulting in a smaller H ii compared to models without dust.

Photodissociation regions
Models including treatment of PDRs are the most computationally intensive in this paper. With these models we can follow ionised hydrogen, through to the atomic and molecular transitions (and can do the same for, e.g. carbon and CO and many other species). The PDR does not affect the temperature in the ionised gas, but leads to heating of up to around a few hundred Kelvin in the neutral gas out to relatively large distances (potentially many parsecs).
In the derivation of Spitzer and H-I equations, there is a term involving the difference in the pressures in the ionised and neutral gas Pi − Po, where Po is dropped assuming Po << Pi. If the PDR heating raises the temperature enough then this assumption will no longer be appropriate and the expansion will be slower. In Bisbas et al. (2015)  respectively. These expressions need to be evaluated numerically. In each case it is the µiTo/µoTi term (which is generally expected to be small) that is dropped to arrive at the Spitzer or Hosokawa-Inutsuka results. PDR (or dust) heating may change this though. The evolution of the H ii region radius for these PDR-RHD models is given in the bottom right panels of Figures 3 and 4. In both density regimes the effect of the PDR is to marginally slow the H ii region expansion. The PDR heated gas close to the ionisation front is ∼ 300 K in both density regimes. Substituting this temperature and an ionised gas temperature of 8200 K into equation 11 and comparing with equation 2 gives an expected difference at 500 kyr of about 1 per cent in the H ii region extent. This is in agreement with the difference found in our simulations. Given the small difference from downstream heating in the PDR model, we also expect the similar heating from dust (see Figure 2) to have a negligible effect.

Consequences for RHD modelling
We have considered the effect of many different microphysical processes that are not normally included in RHD simulations because they make the calculation much more expensive (for the models including PDR treatment, by an order of magnitude or more). Studying all of these processes at once in a systematic way has the added value that they will not be studied across multiple papers and we can immediately see which processes are important and in what ways current simplified models are deficient.

The state of simplified models
Compared to the simplified models in this paper, every other process (except lowering the gas metallicity) results in a smaller H ii region. It is therefore possible that RHD simulations of more complex geometrical systems (such as turbulent star forming regions) using simplified microphysics are overestimating the power of the radiation field. In the context of, for example, radiative feedback and triggered star formation, this could mean that the ability of radiative feedback to compress clouds and trigger star formation might be reduced. Conversely its ability to disperse clouds, halting star formation is also reduced.
The metallicity is also particularly important, with approximately a 20 per cent difference in H ii region radius after 500kyr between models with metallicities 0.5 and 2 times our base metallicity. Radiative feedback may therefore play a substantially different role in the earlier universe or LMC compared to somewhere like the Galactic centre. Equation 10, a simple thermal prescription which accounts for the gas metallicity, could be used to investigate this in 3D models.

Hierarchy of processes when modelling D-type expansion of galactic H ii regions
We have shown here that if only the dynamical evolution of the system needs to be modelled (as opposed to chemistry or synthetic observations) there is not necessarily much value in spending time implementing, for example, detailed spectral models or PDR treatment, when the effect on the dynamical evolution of the system is negligible. We have therefore constructed a simple hierarchy of processes that we find to be important for modelling radiative feedback in star forming regions, which we show in Figure 9. Of course the effects will not be negligible in every scenario (for example radiation pressure is important early on in the formation of massive stars). The hierarchy consists of 4 tiers. Tier 1 is the standard basic RHD (photoionisation + hydrodynamics) code. Tier 2 includes processes that alter the H ii region extent at 500 kyr by around 10 per cent or greater. Tier 2 includes processes that alter the H ii region extent at 500 kyr by of around 1 per cent and tier 4 processes are untested. If improving the microphysics in a photoionisation + hydrodynamics code for non-cosmological applications this should offer a useful reference as to which processes are most important to include. We note that in reality the gas and stellar metallicities are linked. Our hierarchy does not suggest that they are decoupled, rather that if the overall metallicity varies it is the changes in the gas properties rather than the stellar spectrum that affect changes in the dynamical evolution of the system.

Untested features
There are some processes that we identify that might also affect the D-type expansion of galactic H ii regions that we do not treat in this paper. Figure 9. A hierarchy of importance of physical processes for the D-type expansion model considered in this paper. There will be other scenarios where some of these components become significantly more important, e.g. massive star clusters such as the tarantula nebula in the LMC where radiation pressure plays an important role. Note that the percentage difference quoted is only after 500 kyr of H ii region evolution.

X-rays and heavier metals
Treatment of X-ray photons and a wider range of metal species might affect the ionisation and temperature structure in the H ii region. Such treatment has been implemented in the Monte Carlo photoionisation code moccasin (Ercolano et al. 2008) but is not currently available in torus.
Inclusion of X-rays may lead to increased temperatures in the H ii region since atoms/ions with higher ionisation potentials will be ionised. This would lead to an increase in the expansion rate as we found with the lower gas metallicity models.

Non-equilibrium photoionisation
Including non-equilibrium photoionisaiton is not expected to modify the results of the calculations in this paper. It would allow us to capture the R-type expansion of the ionised gas that occurs prior to the D-type expansion that we study here, however R-type expansion occurs rapidly and with little disruption to the density field. Although we do not expect this assumption to affect the resulting density field the results of  show that it will be important to test the effects of non-equilibrium photoionisation on shadowed regions.

SUMMARY AND CONCLUSIONS
We have used the Monte Carlo radiation transport and hydrodynamics code torus and its PDR counterpart torus-3dpdr to study the effects of different microphysics on the D-type expansion of H ii regions. We have compared analytic solutions with the simplest RHD models, calculations including multiple species, the diffuse field, polychromatic radiation, detailed spectral models, different gas metallicity, radiation pressure, dust and PDR's. Each of these processes has been treated directly rather than using a simple parameterisation. We draw the following main conclusions from this work: 1) The critical factor that affects the H ii region dynamics when different microphysics is treated is the gas temperature, which affects the Strömgren radius as rs ∝ T 0.27 and the later time expansion rate as rI ∝ T 2/7 . The Strömgren radius directly affects the expansion rate at times up until 7 c I t 4 rs >> 1.
2) Compared to simplified (hydrogen only, monochromatic, on the spot approximation) models, cooler average gas temperatures resulting from full Monte Carlo photoionisation can change the extent of the H ii region after 500 kyr by about 10 per cent (see conclusion 1). This implies that simplified RHD models are overestimating the power of radiative feedback, be it to induce star formation or disperse gas in star forming regions.
3) Varying the gas metallicity changes the gas amount of metal line cooling and hence the gas temperature. Changing the metallicity by a factor of 4 can also lead to a 10-20 per cent difference in the H ii region extent after 500 kyr. The Spitzer D-type expansion solution matches our metallicity dependent models well if it uses the average gas temperatures from our simulation. We fit our models to produce a simplified thermal prescription that can be used in 3D models to incorporate heating effects of multi-species gas at varying metallicity. This is at negligible additional computational cost. It will also allow us to test how deficient simplified schemes are in more complex 3D applications.
3) Dust absorption reduces the ionising photon budget, meaning that the onset of D-type expansion occurs at smaller radii (since rs ∝ N 1/3 ly ). The H ii region temperature is similar, so the later time expansion rate (when the Strömgren radius term becomes less important) is similar with or without dust. 4) Using detailed stellar spectral models does not affect the simulation much (of order 1 per cent difference in the ionised gas extent at 500 kyr), since the ionising flux remains approximately constant and the effect of any harder radiation is just to smear out the ionisation front. Radiation pressure is also found to be dominated by the effects of photoionisation, in agreement with previous studies such as Sales et al. (2014). PDR heating of the downstream medium only marginally affects the H ii region expansion (by ∼ 1 per cent at 500 kyr), in agreement with the difference expected from the analytic equations. 6) We develop a hierarchy of processes for modelling the D-type expansion of H ii regions, where additional physics beyond the most simple model is graded based upon its effect upon the expansion. This offers guidance for where to focus in future development of numerical models.

ACKNOWLEDGMENTS
TJHaw is funded by the STFC consolidated grant ST/K000985/1. TJHarries acknowledges funding via an STFC consolidated grant ST/J001627/1. The work of TGB was funded by STFC grant ST/J001511/1.
The calculations presented here were performed using the Darwin Data Analytic system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a BIS National Einfrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure. This work was also supported by initial exploratory calculations performed using the DiRAC HPC facility machines Zen (Exeter) and Complexity (Leicester).
We would like to acknowledge the Nordita program on Photo-Evaporation in Astrophysical Systems (June 2013) where part of the work for this paper was carried out. This work also benefitted from discussions held at the Starbench (Exeter, 2013) and Starbench-2 (Bönn, 2014) code comparison project meetings. We also thank Barbara Ercolano for her assistance in comparing torus with mocassin in the dust decoupled HII40 Lexington benchmark.