The SILCC project

We present magnetohydrodynamic (MHD) simulations of the star-forming multiphase interstellar medium (ISM) in stratified galactic patches with gas surface densities $\Sigma_\mathrm{gas} = 10, 30, 50, 100\,\mathrm{M_\odot\,pc^{-2}}$. The SILCC project simulation framework accounts for non-equilibrium thermal and chemical processes in the warm and cold ISM. The sink-based star formation and feedback model includes stellar winds, hydrogen-ionising UV radiation, core-collapse supernovae, and cosmic ray (CR) injection and diffusion. The simulations follow the observed relation between $\Sigma_\mathrm{gas}$ and the star formation rate surface density, $\Sigma_\mathrm{SFR}$. CRs qualitatively change the outflow phase structure. Without CRs, the outflows transition from a two-phase (warm and hot at 1 kpc) to a single-phase (hot at 2 kpc) structure. With CRs, the outflow always has three phases (cold, warm, and hot), dominated in mass by the warm phase. The impact of CRs on mass loading decreases for higher $\Sigma_\mathrm{gas}$ and the mass loading factors of the CR-supported outflows are of order unity independent of $\Sigma_\mathrm{SFR}$. Similar to observations, vertical velocity dispersions of the warm ionised medium (WIM) and the cold neutral medium (CNM) correlate with the star formation rate as $\sigma_\mathrm{z} \propto \Sigma^{\alpha}_\mathrm{SFR}$, with $\alpha \sim 0.20$. In the absence of stellar feedback, we find no correlation. The velocity dispersion of the WIM is a factor $\sim 2.2$ higher than that of the CNM, in agreement with local observations. For $\Sigma_\mathrm{SFR} \geq 1.5 \times 10^{-2}\,\mathrm{M_\odot\,yr^{-1}\,kpc^{-2}}$ the WIM motions become supersonic.

This paper was published in the Monthly Notices of the Royal Astronomical Society (April 2023), Volume 522, Issue 2, p. 1843-1862. A copy of the preprint is also available on the astro-ph server (arXiv:2211.15419).

We present simulations of the multiphase interstellar medium (ISM) at solar neighbourhood conditions including thermal and non-thermal ISM processes, star cluster formation, and feedback from massive stars: stellar winds, hydrogen ionizing radiation computed with the novel TREERAY radiative transfer method, supernovae (SN), and the injection of cosmic rays (CR). N-body dynamics is computed with a 4th-order Hermite integrator. We systematically investigate the impact of stellar feedback on the self-gravitating ISM with magnetic fields, CR advection and diffusion, and non-equilibrium chemical evolution. SN-only feedback results in strongly clustered star formation with very high star cluster masses, a bi-modal distribution of the ambient SN densities, and low volume-filling factors (VFF) of warm gas, typically inconsistent with local conditions. Early radiative feedback prevents an initial starburst, reduces star cluster masses and outflow rates. Furthermore, star formation rate surface densities of $\Sigma_\mathrm{SFR} = 1.4 - 5.9 \times 10^{-3} \,\mathrm{M}_\odot\,\mathrm{yr}^{-1}\,\mathrm{kpc}^{-2}$, VFF$_\mathrm{warm} = 60 - 80$ per cent as well as thermal, kinetic, magnetic, and cosmic ray energy densities of the model including all feedback mechanisms agree well with observational constraints. On the short, 100 Myr, time-scales investigated here, CRs only have a moderate impact on star formation and the multiphase gas structure and result in cooler outflows, if present. Our models indicate that at low gas surface densities SN-only feedback only captures some characteristics of the star-forming ISM and outflows/inflows relevant for regulating star formation. Instead, star formation is regulated on star cluster scales by radiation and winds from massive stars in clusters, whose peak masses agree with solar neighbourhood estimates.

This paper was published in the Monthly Notices of the Royal Astronomical Society (June 2021), Volume 504, Issue 1, p. 1039-1061. A copy of the preprint is also available on the astro-ph server (arXiv:2103.14128).

Magnetic fields are ubiquitously observed in the interstellar medium (ISM) of present-day star-forming galaxies with dynamically relevant energy densities. Using three-dimensional magnetohydrodynamic (MHD) simulations of the supernova-driven ISM in the flux-freezing approximation (ideal MHD), we investigate the impact of the magnetic field on the chemical and dynamical evolution of the gas, fragmentation, and the formation of molecular clouds. We follow the chemistry with a network of six species ($H^+$, $H$, $H_2$, $C^+$, $CO$, and free electrons) including local shielding effects. We find that magnetic fields thicken the disc by a factor of a few to a scale height of $\sim 100$ pc⁠, delay the formation of dense (and molecular) gas by $\sim 25$ Myr⁠, and result in differently shaped gas structures. The magnetized gas fragments into fewer clumps, which are initially at subcritical mass-to-flux ratios, $M/\Phi \approx 0.3(M/Φ)_{crit}$, and accrete gas preferentially parallel to the magnetic field lines until supercritial mass-to-flux ratios of up to the order of 10 are reached. The accretion rates onto molecular clouds scale with $\dot{M} \propto M^{1.5}$⁠. The median of the intercloud velocity dispersion is $\sim 2−5$ $kms^{−1}$ and lower than the internal velocity dispersion in the clouds (⁠$\sim 3−7$ $kms^{−1}$⁠). However, individual cloud–cloud collisions occur at speeds of a few $10$ $kms^{−1}$⁠.

This paper was published in the Monthly Notices of the Royal Astronomical Society (November 1, 2018), Volume 480, Issue 3, p. 3511-3540. A copy of the preprint is also available on the astro-ph server (arXiv:1808.05222).

Dust grains are an important component of the interstellar medium (ISM) of galaxies. We present the first direct measurement of the residence times of interstellar dust in the different ISM phases, and of the transition rates between these phases, in realistic hydrodynamical simulations of the multi-phase ISM. Our simulations include a time-dependent chemical network that follows the abundances of H${}^+$, H, H${}_2$, C${}^+$ and CO and take into account self-shielding by gas and dust using a tree-based radiation transfer method. Supernova explosions are injected either at random locations, at density peaks, or as a mixture of the two. For each simulation, we investigate how matter circulates between the ISM phases and find more sizeable transitions than considered in simple mass exchange schemes in the literature. The derived residence times in the ISM phases are characterised by broad distributions, in particular for the molecular, warm and hot medium. The most realistic simulations with random and mixed driving have median residence times in the molecular, cold, warm and hot phase around 17, 7, 44 and 1 Myr, respectively. The transition rates measured in the random driving run are in good agreement with observations of Ti gas-phase depletion in the warm and cold phases in a simple depletion model, although the depletion in the molecular phase is under-predicted. ISM phase definitions based on chemical abundance rather than temperature cuts are physically more meaningful, but lead to significantly different transition rates and residence times because there is no direct correspondence between the two definitions.

This paper has been published in the Monthly Notices of the Royal Astronomical Society (June 1, 2017), Volume 467, Issue 4, p.4322-4342. A copy of the preprint is also available on the astro-ph server (arXiv:1610.06579).

We present three-dimensional radiation-hydrodynamical simulations of the impact of stellar winds, photoelectric heating, photodissociating and photoionising radiation, and supernovae on the chemical composition and star formation in a stratified disc model. This is followed with a sink-based model for star clusters with populations of individual massive stars. Stellar winds and ionising radiation regulate the star formation rate at a factor of ~10 below the simulation with only supernova feedback due to their immediate impact on the ambient interstellar medium after star formation. Ionising radiation (with winds and supernovae) significantly reduces the ambient densities for most supernova explosions to rho < 10^-25 g cm^-3, compared to 10^-23 g cm^-3 for the model with only winds and supernovae. Radiation from massive stars reduces the amount of molecular hydrogen and increases the neutral hydrogen mass and volume filling fraction. Only this model results in a molecular gas depletion time scale of 2 Gyr and shows the best agreement with observations. In the radiative models, the Halpha emission is dominated by radiative recombination as opposed to collisional excitation (the dominant emission in non-radiative models), which only contributes ~1-10 % to the total Halpha emission. Individual massive stars ($M>= 30 M_\odot$) with short lifetimes are responsible for significant fluctuations in the Halpha luminosities. The corresponding inferred star formation rates can underestimate the true instantaneous star formation rate by factors of ~10.

This paper has been published in the Monthly Notices of the Royal Astronomical Society (December 10, 2016), Volume 466, Issue 3, p. 3293-3308. A copy of the preprint is also available on the astro-ph server (arXiv:1610.06569).

We study the impact of stellar winds and supernovae on the multi-phase interstellar medium using three-dimensional hydrodynamical simulations carried out with FLASH. The selected galactic disc region has a size of (500 pc)$^2$ x $\pm$ 5 kpc and a gas surface density of 10 M$_{\odot}$/pc$^2$. The simulations include an external stellar potential and gas self-gravity, radiative cooling and diffuse heating, sink particles representing star clusters, stellar winds from these clusters which combine the winds from indi- vidual massive stars by following their evolution tracks, and subsequent supernova explosions. Dust and gas (self-)shielding is followed to compute the chemical state of the gas with a chemical network. We find that stellar winds can regulate star (cluster) formation. Since the winds suppress the accretion of fresh gas soon after the cluster has formed, they lead to clusters which have lower average masses (10$^2$ - 10$^{4.3}$ M$_{\odot}$) and form on shorter timescales (10$^{-3}$ - 10 Myr). In particular we find an anti-correlation of cluster mass and accretion time scale. Without winds the star clusters easily grow to larger masses for ~5 Myr until the first supernova explodes. Overall the most massive stars provide the most wind energy input, while objects beginning their evolution as B-type stars contribute most of the supernova energy input. A significant outflow from the disk (mass loading $\gtrsim$ 1 at 1 kpc) can be launched by thermal gas pressure if more than 50% of the volume near the disc mid-plane can be heated to T > 3x10$^5$ K. Stellar winds alone cannot create a hot volume-filling phase. The models which are in best agreement with observed star formation rates drive either no outflows or weak outflows.

This paper has been published in the Monthly Notices of the Royal Astronomical Society (December 12, 2016), Volume 466, Issue 2, p. 1903-1924. A copy of the preprint is also available on the astro-ph server (arXiv:1606.05346).

The SILCC project (SImulating the Life-Cycle of molecular Clouds) aims at a more self-consistent understanding of the interstellar medium (ISM) on small scales and its link to galaxy evolution. We present three-dimensional (magneto)hydrodynamic simulations of the ISM in a vertically stratified box including self-gravity, an external potential due to the stellar component of the galactic disc, and stellar feedback in the form of an interstellar radiation field and supernovae (SNe). The cooling of the gas is based on a chemical network that follows the abundances of ${\rm H}^{+}$, ${\rm H}$, ${\rm H}_2$, ${\rm C}^{+}$, and ${\rm CO}$ and takes shielding into account consistently. We vary the SN feedback by comparing different SN rates, clustering and different positioning, in particular SNe in density peaks and at random positions, which has a major impact on the dynamics. Only for random SN positions the energy is injected in sufficiently low-density environments to reduce energy losses and enhance the effective kinetic coupling of the SNe with the gas. This leads to more realistic velocity dispersions ($\sigma_{HI} $$~\sim~$$ 0.8\sigma_{(300 - 8000{\rm~K})} $$~\sim~$$ 10-20{\rm~km/s}$, $\sigma_{H\alpha} $$~\sim~$$ 0.6\sigma_{(8000 - 3 \cdot 10^5 {\rm~K})} $$~\sim~$$ 20-30{\rm~km/s}$), and strong outflows with mass loading factors of up to 10 even for solar neighbourhood conditions. Clustered SNe abet the onset of outflows compared to individual SNe but do not influence the net outflow rate. The outflows do not contain any molecular gas and are mainly composed of atomic hydrogen. The bulk of the outflowing mass is dense ($\rho $$~\sim~$$ 10^{-25} - 10^{-24} {\rm~g/cc}$) and slow ($v $$~\sim~$$ 20-40 {\rm~km/s}$) but there is a high-velocity tail of up to $v $$~\sim~$$ 500{\rm~km/s}$ with $\rho $$~\sim~$$ 10^{-28} - 10^{-27} {\rm~g/cc}$.

This paper has been published in the Monthly Notices of the Royal Astronomical Society (March 11, 2016) 456 (4): 3432-3455. A copy of the preprint is also available on the astro-ph server (arXiv:1508.06646).

We present a hydrodynamical simulation of the turbulent, magnetized, supernova (SN)-driven interstellar medium (ISM) in a stratified box that dynamically couples the injection and evolution of cosmic rays (CRs) and a self-consistent evolution of the chemical composition. CRs are treated as a relativistic fluid in the advection-diffusion approximation. The thermodynamic evolution of the gas is computed using a chemical network that follows the abundances of H${\rm H}^{+}$, ${\rm H}$, ${\rm H}_2$, ${\rm CO}$, and ${\rm C}^{+}$, and free electrons and includes (self-)shielding of the gas and dust. We find that CRs perceptibly thicken the disk with the heights of 90% (70%) enclosed mass reaching ~1.5 kpc (~0.2 kpc). The simulations indicate that CRs alone can launch and sustain strong outflows of atomic and ionized gas with mass loading factors of order unity, even in solar neighborhood conditions and with a CR energy injection per SN of 10^50 erg, 10% of the fiducial thermal energy of an SN. The CR-driven outflows have moderate launching velocities close to the midplane (~100 km/s) and are denser ($\rho$~1e-24 - 1e-26 g/cm^3), smoother, and colder than the (thermal) SN-driven winds. The simulations support the importance of CRs for setting the vertical structure of the disk as well as the driving of winds.

This paper has been published in the The Astrophysical Journal Letters, Volume 816, Issue 2. A copy of the preprint is also available on the astro-ph server (arXiv:1509.07247).

The SILCC (SImulating the Life-Cycle of molecular Clouds) project aims to self-consistently understand the small-scale structure of the interstellar medium (ISM) and its link to galaxy evolution. We simulate the evolution of the multiphase ISM in a $(500{\rm~pc})^2 \times \pm 5 {\rm~kpc}$ region of a galactic disc, with a gas surface density of $\Sigma _{_{\rm GAS}} = 10 \;{\rm M}_{\odot }\,{\rm pc}^{-2}$. The $\mathtt{FLASH}$ 4 simulations include an external potential, self-gravity, magnetic fields, heating and radiative cooling, time-dependent chemistry of ${\rm H}_2$ and ${\rm CO}$ considering (self-) shielding, and supernova (SN) feedback but omit shear due to galactic rotation. We explore SN explosions at different rates in high-density regions (peak), in random locations with a Gaussian distribution in the vertical direction (random), in a combination of both (mixed), or clustered in space and time (clus/clus2). Only models with self-gravity and a significant fraction of SNe that explode in low-density gas are in agreement with observations. Without self-gravity and in models with peak driving the formation of ${\rm H}_2$ is strongly suppressed. For decreasing SN rates, the ${\rm H}_2$ mass fraction increases significantly from $< 10$ per cent for high SN rates, i.e. 0.5 dex above Kennicutt–Schmidt, to 70–85 per cent for low SN rates, i.e. 0.5 dex below KS. For an intermediate SN rate, clustered driving results in slightly more ${\rm H}_2$ than random driving due to the more coherent compression of the gas in larger bubbles. Magnetic fields have little impact on the final disc structure but affect the dense gas ($n \gtrsim 10{\rm~cm}^{-3}$) and delay ${\rm H}_2$ formation. Most of the volume is filled with hot gas ($\sim 80$ per cent within $\pm 150{\rm~pc}$). For all but peak driving a vertically expanding warm component of atomic hydrogen indicates a fountain flow. We highlight that individual chemical species populate different ISM phases and cannot be accurately modelled with temperature-/density-based phase cut-offs.

This paper has been published in the Monthly Notices of the Royal Astronomical Society (November 21, 2015) 454 (1): 246-276. A copy of the preprint is also available on the astro-ph server (arXiv:1412.2749).

We present a study of the cold atomic hydrogen (HI) content of molecular clouds simulated within the SILCC-Zoom project for solar neighbourhood conditions. We produce synthetic observations of HI at 21 cm including HI self-absorption (HISA) and observational effects. We find that HI column densities, NHI , of ≳ 1022 cm-2 are frequently reached in molecular clouds with HI temperatures as low as ~10 K. Hence, HISA observations assuming a fixed HI temperature tend to underestimate the amount of cold HI in molecular clouds by a factor of 3 - 10 and produce an artificial upper limit of NHI around 1021 cm-2. We thus argue that the cold HI mass in molecular clouds could be a factor of a few higher than previously estimated. Also NHI -PDFs obtained from HISA observations might be subject to observational biases and should be considered with caution. The underestimation of cold HI in HISA observations is due to both the large HI temperature variations and the effect of noise in regions of high optical depth. We find optical depths of cold HI around 1 - 10 making optical depth corrections essential. We show that the high HI column densities (≳ 1022 cm-2) can in parts be attributed to the occurrence of up to 10 individual HI-H2 transitions along the line of sight. This is also reflected in the spectra, necessitating Gaussian decomposition algorithms for their in-depth analysis. However, also for a single HI-H2 transition, NHI frequently exceeds 1021 cm-2, challenging 1-dimensional, semi-analytical models. This is due to non-equilibrium chemistry effects and the fact that HI-H2 transition regions usually do not possess a 1-dimensional geometry. Finally, we show that the HI gas is moderately supersonic with Mach numbers of a few. The corresponding non-thermal velocity dispersion can be determined via HISA observations within a factor of ~2.

This paper was published in the Monthly Notices of the Royal Astronomical Society, Advance Access. A copy of the preprint is also available on the astro-ph server (arXiv:2109.10917).

The properties of synthetic CO emission from 3D simulations of forming molecular clouds are studied within the SILCC-Zoom project. Since the time-scales of cloud evolution and molecule formation are comparable, the simulations include a live chemical network. Two sets of simulations with an increasing spatial resolution (dx = 3.9 pc to dx = 0.06 pc) are used to investigate the convergence of the synthetic CO emission, which is computed by post-processing the simulation data with the RADMC-3D radiative transfer code. To determine the excitation conditions, it is necessary to include atomic hydrogen and helium alongside H2, which increases the resulting CO emission by ~7-26 per cent. Combining the brightness temperature of 12CO and 13CO, we compare different methods to estimate the excitation temperature, the optical depth of the CO line and hence, the CO column density. An intensity-weighted average excitation temperature results in the most accurate estimate of the total CO mass. When the pixel-based excitation temperature is used to calculate the CO mass, it is over-/underestimated at low/high CO column densities where the assumption that 12CO is optically thick while 13CO is optically thin is not valid. Further, in order to obtain a converged total CO luminosity and hence ⟨XCO⟩ factor, the 3D simulation must have dx ≲ 0.1 pc. The ⟨XCO⟩ evolves over time and differs for the two clouds; yet pronounced differences with numerical resolution are found. Since high column density regions with a visual extinction larger than 3 mag are not resolved for dx ≳ 1 pc, in this case the H2 mass and CO luminosity both differ significantly from the higher resolution results and the local XCO is subject to strong noise. Our calculations suggest that synthetic CO emission maps are only converged for simulations with dx ≲ 0.1 pc.

This paper was published in the Monthly Notices of the Royal Astronomical Society, Volume 510, Issue 1, pp.753-773. A copy of the preprint is also available on the astro-ph server (arXiv:2102.00778).

We present synthetic dust polarization maps of simulated molecular clouds with the goal to systematically explore the origin of the relative orientation of the magnetic field ( B ) with respect to the cloud sub-structure identified in density (n; 3D) and column density (N; 2D). The polarization maps are generated with the radiative transfer code POLARIS, which includes self-consistently calculated efficiencies for radiative torque alignment. The molecular clouds are formed in two sets of 3D magnetohydrodynamical simulations: (i) in colliding flows (CF), and (ii) in the SILCC-Zoom simulations. In 3D, for the CF simulations with an initial field strength below ∼5 μG, B is oriented either parallel or randomly with respect to the n-structures. For CF runs with stronger initial fields as well as all SILCC-Zoom simulations, which have an initial field strength of 3 μG, a flip from parallel to perpendicular orientation occurs at high densities of ntrans≃ 102-103 cm-3. We suggest that this flip happens if the cloud's mass-to-flux ratio, μ, is close to or below the critical value of 1. This corresponds to a field strength around 3-5 μG, close to the Galactic average. In 2D, we use the method of Projected Rayleigh Statistics (PRS) to study the relative orientation of B . If present, the flip in orientation occurs in the projected maps at Ntrans≃ 1021 - 21.5 cm-2. This value is similar to the observed transition value from sub- to supercritical magnetic fields in the interstellar medium. However, projection effects can strongly reduce the predictive power of the PRS method: Depending on the considered cloud or line-of-sight, the projected maps of the SILCC-Zoom simulations do not always show the flip, although it is expected given the 3D morphology. Such projection effects can explain the variety of recently observed field configurations, in particular within a single cloud. Finally, we do not find a correlation between the observed orientation of B and the N-PDF.

This paper was published in the Monthly Notices of the Royal Astronomical Society, Volume 497, Issue 4, pp.4196-4212. A copy of the preprint is also available on the astro-ph server (arXiv:2003.00017).

As part of the SILCC-Zoom project, we present our first sub-parsec resolution radiation-hydrodynamic simulations of two molecular clouds self-consistently forming from a turbulent, multiphase ISM. The clouds have similar initial masses of few $10^4 M_{\odot}$, escape velocities of $5 km s^{-1}$, and a similar initial energy budget. We follow the formation of star clusters with a sink-based model and the impact of radiation from individual massive stars with the tree-based radiation transfer module TREERAY. Photoionizing radiation is coupled to a chemical network to follow gas heating, cooling, and molecule formation and dissociation. For the first 3 Myr of cloud evolution, we find that the overall star formation efficiency is considerably reduced by a factor of $4$ to global cloud values of $\leq 10$ per cent as the mass accretion of sinks that host massive stars is terminated after $\leq 1$ Myr. Despite the low efficiency, star formation is triggered across the clouds. Therefore, a much larger region of the cloud is affected by radiation and the clouds begin to disperse. The time-scale on which the clouds are dispersed sensitively depends on the cloud sub-structure and in particular on the amount of gas at high visual extinction. The damage of radiation done to the highly shielded cloud ($MC_1$) is delayed. We also show that the radiation input can sustain the thermal and kinetic energy of the clouds at a constant level. Our results strongly support the importance of ionizing radiation from massive stars for explaining the low-observed star formation efficiency of molecular clouds.

This paper was published in the Monthly Notices of the Royal Astronomical Society (December 19, 2019), Volume 492, Issue 1, p. 1465-1483. A copy of the preprint is also available on the astro-ph server (arXiv:1906.01015).

As part of the SILCC-Zoom project, we present our first sub-parsec resolution radiation-hydrodynamic simulations of two molecular clouds self-consistently forming from a turbulent, multiphase ISM. The clouds have similar initial masses of few $10^4 M_{\odot}$, escape velocities of $5 km s^{-1}$, and a similar initial energy budget. We follow the formation of star clusters with a sink-based model and the impact of radiation from individual massive stars with the tree-based radiation transfer module TREERAY. Photoionizing radiation is coupled to a chemical network to follow gas heating, cooling, and molecule formation and dissociation. For the first 3 Myr of cloud evolution, we find that the overall star formation efficiency is considerably reduced by a factor of $4$ to global cloud values of $\leq 10$ per cent as the mass accretion of sinks that host massive stars is terminated after $\leq 1$ Myr. Despite the low efficiency, star formation is triggered across the clouds. Therefore, a much larger region of the cloud is affected by radiation and the clouds begin to disperse. The time-scale on which the clouds are dispersed sensitively depends on the cloud sub-structure and in particular on the amount of gas at high visual extinction. The damage of radiation done to the highly shielded cloud ($MC_1$) is delayed. We also show that the radiation input can sustain the thermal and kinetic energy of the clouds at a constant level. Our results strongly support the importance of ionizing radiation from massive stars for explaining the low-observed star formation efficiency of molecular clouds.

This paper was published in the Monthly Notices of the Royal Astronomical Society (January 21, 2019), Volume 482, Issue 3, p. 4062-4083. A copy of the preprint is also available on the astro-ph server (arXiv:1810.08210).

We present synthetic dust polarisation maps of 3D magneto-hydrodynamical simulations of molecular clouds embedded in their galactic environment performed within the SILCC-Zoom project. The radiative transfer is carried out with POLARIS for wavelengths from 70 $\mu$m to 3 mm at a resolution of 0.12 pc, and includes self-consistently calculated alignment efficiencies for radiative torque alignment. We explore the reason of the observed depolarisation in the center of molecular clouds: We find that dust grains remain well aligned even at high densities (n > 103 cm$^{−3}$) and visual extinctions (A$_V$ > 1). The depolarisation is rather caused by strong variations of the magnetic field direction along the LOS due to turbulent motions. The observed magnetic field structure thus resembles best the mass-weighted, line-of-sight averaged field structure. Furthermore, it differs by only a few 1$^{\circ}$ for different wavelengths and is little affected by the spatial resolution of the synthetic observations. Noise effects can be reduced by convolving the image. Doing so, for $\lambda \geq 160$ $\mu$m the observed magnetic field traces reliably the underlying field in regions with intensities I $\geq$ 3 times the noise level and column densities above 1 M$_{sun}$ pc$^{−2}$. Here, typical deviations are $\leq$ 10$^{\circ}$. The observed structure is less reliable in regions with low polarisation degrees and possibly in regions with large column density gradients. Finally, we show that the simplified and widely used method by Wardle & Konigl (1990) without self-consistent dust alignment efficiencies can provide a good representation of the observable field structure with deviations below 5$^{\circ}$.

This paper was published in the Monthly Notices of the Royal Astronomical Society (January 11, 2019), Volume 482, Issue 2, p. 2697-2716. A copy of the preprint is also available on the astro-ph server (arXiv:1804.10157).

The C+ ion is an important coolant of interstellar gas, and so the [C II] fine structure line is frequently observed in the interstellar medium. However, the physical and chemical properties of the [C II]-emitting gas are still unclear. We carry out non-LTE (local thermal equilibrium) radiative transfer simulations with RADMC-3D to study the [C II] line emission from a young, turbulent molecular cloud before the onset of star formation, using data from the SILCC-Zoom project. The [C II] emission is optically thick over 40 per cent of the observable area with I[C II] > 0.5 K km s-1. To determine the physical properties of the [C II] emitting gas, we treat the [C II] emission as optically thin. We find that the [C II] emission originates primarily from cold, moderate density gas (40 ≲ T ≲ 65 K and 50 ≲ n ≲ 440 cm-3), composed mainly of atomic hydrogen and with an effective visual extinction between ∼0.50 and ∼0.91. Gas dominated by molecular hydrogen contributes only ≲20 per cent of the total [C II] line emission. Thus, [C II] is not a good tracer for CO-dark H2 at this early phase in the cloud's lifetime. We also find that the total gas, H and C+ column densities are all correlated with the integrated [C II] line emission, with power law slopes ranging from 0.5 to 0.7. Further, the median ratio between the total column density and the [C II] line emission is YC II ≈ 1.1 × 1021 cm-2 (K km s-1)-1, and Y_{C II}} scales with I_{[C II]}^{-0.3}. We expect Y_{C II} to change in environments with a lower or higher radiation field than simulated here.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 481, Issue 4, p.4277-4299 (2018). A copy of the preprint is also available on the astro-ph server (arXiv:1809.10696).

We present high-resolution ($\sim$0.1 pc), hydrodynamical and magnetohydrodynamical simulations to investigate whether the observed level of molecular cloud (MC) turbulence can be generated and maintained by external supernova (SN) explosions. The MCs are formed self-consistently within their large-scale galactic environment following the non-equilibrium formation of H2 and CO, including (self-) shielding and important heating and cooling processes. The MCs inherit their initial level of turbulence from the diffuse ISM, where turbulence is injected by SN explosions. However, by systematically exploring the effect of individual SNe going off outside the clouds, we show that at later stages the importance of SN-driven turbulence is decreased significantly. This holds for different MC masses as well as for MCs with and without magnetic fields. The SN impact also decreases rapidly with larger distances. Nearby SNe (d $\sim$ 25 pc) boost the turbulent velocity dispersions of the MC by up to 70% (up to a few km s$^{‑1}$). For d $\geq$ 50 pc, however, their impact decreases fast with increasing d and is almost negligible. For all probed distances the gain in velocity dispersion decays rapidly within a few 100 kyr. This is significantly shorter than the average timescale for an MC to be hit by a nearby SN under solar neighborhood conditions ($\sim$2 Myr). Hence, at these conditions SNe are not able to sustain the observed level of MC turbulence. However, in environments with high gas surface densities and SN rates, like the Central Molecular Zone, observed elevated MC dispersions could be triggered by external SNe.

This paper has been published in the The Astrophysical Journal, Volume 855, Issue 2, article id. 81, 9 pp. (2018). A copy of the preprint is also available on the astro-ph server (arXiv:1802.00973).

We present 3D zoom-in simulations of the formation of two molecular clouds out of the galactic interstellar medium. We model the clouds – identified from the SILCC simulations – with a resolution of up to 0.06 pc using adaptive mesh refinement in combination with a chemical network to follow heating, cooling and the formation of $\rm~H_2$ and $\rm~CO$ including (self-) shielding. The two clouds are assembled within a few million years with mass growth rates of up to $\sim 10^{−2}$ ${\rm M}_{\odot} \rm~yr^{−1}$ and final masses of $\sim 50000 \rm~M_{\odot}$. A spatial resolution of $\leqslant 0.1 \rm~pc$ is required for convergence with respect to the mass, velocity dispersion and chemical abundances of the clouds, although these properties also depend on the cloud definition such as based on density thresholds, $\rm~H_2$ or $\rm~CO$ mass fraction. To avoid grid artefacts, the progressive increase of resolution has to occur within the free-fall time of the densest structures ($1–1.5 \rm~Myr$) and $\geqslant 200$ time-steps should be spent on each refinement level before the resolution is progressively increased further. This avoids the formation of spurious, large-scale, rotating clumps from unresolved turbulent flows. While $\rm~CO$ is a good tracer for the evolution of dense gas with number densities $\rm{n} \geqslant 300 \rm~cm^{−3}$, $\rm~H_2$ is also found for $\rm{n} \leqslant 30 \rm~cm^{−3}$ due to turbulent mixing and becomes dominant at column densities around $30–50 \rm{M}_{\odot} \rm{pc}^{−2}$. The $\rm~CO$-to-$\rm~H_2$ ratio steadily increases within the first $2 \rm~Myr$, whereas $\rm{X_{CO}} \simeq 1–4 \times 10^{20} \rm~cm^{−2} (K\ km\ s^{−1})^{−1}$ is approximately constant since the CO(1−0) line quickly becomes optically thick.

This paper was published in the Monthly Notices of the Royal Astronomical Society (December 21, 2017), Volume 472, Issue 4, p. 4797-4818. A copy of the preprint is also available on the astro-ph server (arXiv:1704.06487).

The degree to which the formation and evolution of clouds and filaments in the interstellar medium is regulated by magnetic fields remains an open question. Yet the fundamental properties of the fields (strength and 3D morphology) are not readily observable. We investigate the potential for recovering magnetic field information from dust polarization, the Zeeman effect, and the Faraday rotation measure (RM) in a SILCC-Zoom magnetohydrodynamic (MHD) filament simulation. The object is analysed at the onset of star formation and it is characterized by a line-mass of about $(M/L) \sim 63 \mathrm{M_\odot} pc^{-1}$ out to a radius of 1 pc and a kinked 3D magnetic field morphology. We generate synthetic observations via POLARIS radiative transfer (RT) post-processing and compare with an analytical model of helical or kinked field morphology to help interpreting the inferred observational signatures. We show that the tracer signals originate close to the filament spine. We find regions along the filament where the angular dependence with the line of sight (LOS) is the dominant factor and dust polarization may trace the underlying kinked magnetic field morphology. We also find that reversals in the recovered magnetic field direction are not unambiguously associated to any particular morphology. Other physical parameters, such as density or temperature, are relevant and sometimes dominant compared to the magnetic field structure in modulating the observed signal. We demonstrate that the Zeeman effect and the RM recover the line-of-sight magnetic field strength to within a factor 2.1-3.4. We conclude that the magnetic field morphology may not be unambiguously determined in low-mass systems by observations of dust polarization, Zeeman effect, or RM, whereas the field strengths can be reliably recovered.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 500, Issue 1, pp.153-176. A copy of the preprint is also available on the astro-ph server (arXiv:2009.04201).

We take advantage of a set of molecular cloud simulations to demonstrate a possibility to uncover statistical properties of the gas density and velocity fields using reflected emission of a short (with duration much less than the cloud's light-crossing time) X-ray flare. Such a situation is relevant for the Central Molecular Zone (CMZ) of our Galaxy where several clouds get illuminated by an ∼110 yr-old flare from the supermassive black hole Sgr A* . Due to shortness of the flare (Δt ≲ 1.6 yr), only a thin slice (Δz ≲ 0.5 pc) of the molecular gas contributes to the X-ray reflection signal at any given moment, and its surface brightness effectively probes the local gas density. This allows reconstructing the density probability distribution function over a broad range of scales with virtually no influence of attenuation, chemo-dynamical biases, and projection effects. Such a measurement is key to understanding the structure and star formation potential of the clouds evolving under extreme conditions in the CMZ. For cloud parameters similar to the currently brightest in X-ray reflection molecular complex Sgr A, the sensitivity level of the best available data is sufficient only for marginal distinction between solenoidal and compressive forcing of turbulence. Future-generation X-ray observatories with large effective area and high spectral resolution will dramatically improve on that by minimizing systematic uncertainties due to contaminating signals. Furthermore, measurement of the iron fluorescent line centroid with sub-eV accuracy in combination with the data on molecular line emission will allow direct investigation of the gas velocity field.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 495, Issue 1, pp.1414-1432. A copy of the preprint is also available on the astro-ph server (arXiv:1906.11579).

Photoionizing radiation and stellar winds from massive stars deposit energy and momentum into the interstellar medium (ISM). They might disperse the local ISM, change its turbulent multi-phase structure, and even regulate star formation. Ionizing radiation dominates the massive stars' energy output, but the relative effect of winds might change with stellar mass and the properties of the ambient ISM. We present simulations of the interaction of stellar winds and ionizing radiation of 12, 23, and 60 M$_{\odot}$ stars within a cold neutral (CNM, n$_0$ = 100 cm$^{−3}$), warm neutral (WNM, n$_0$ = 1, 10 cm$^{−3}$) or warm ionized (WIM, n$_0$ = 0.1 cm$^{−3}$) medium. The FLASH simulations adopt the novel tree-based radiation transfer algorithm TreeRay. With the On-the-Spot approximation and a temperature-dependent recombination coefficient, it is coupled to a chemical network with radiative heating and cooling. In the homogeneous CNM, the total momentum injection ranges from $1.6 \times 10^4$ to $4 \times 10^5$ M$_{\odot}$ km s$_{−1}$ and is always dominated by the expansion of the ionized HII region. In the WIM, stellar winds dominate ($2 \times 10^2$ to $5 \times 10^3$ M$_{\odot}$ km s$^{−1}$), while the input from radiation is small ($\sim 10^2$ M$_{\odot}$ km s$^{−1}$). The WNM (n$_0$ = 1 cm$^{−3}$) is a transition regime. Energetically, stellar winds couple more efficiently to the ISM ($\sim$ 0.1 percent of wind luminosity) than radiation (< 0.001 percent of ionizing luminosity). For estimating the impact of massive stars, the strongly mass-dependent ratios of wind to ionizing luminosity and the properties of the ambient medium have to be considered.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 478, Issue 4, p.4799-4815 (2018). A copy of the preprint is also available on the astro-ph server (arXiv:1804.10218).

It has been known for more than 30 yr that the distribution of molecular gas in the innermost 300 parsecs of the Milky Way, the Central Molecular Zone, is strongly asymmetric. Indeed, approximately three quarters of molecular emission come from positive longitudes, and only one quarter from negative longitudes. However, despite much theoretical effort, the origin of this asymmetry has remained a mystery. Here, we show that the asymmetry can be neatly explained by unsteady flow of gas in a barred potential. We use high-resolution 3D hydrodynamical simulations coupled to a state-of-the-art chemical network. Despite the initial conditions and the bar potential being point symmetric with respect to the Galactic Centre, asymmetries develop spontaneously due to the combination of a hydrodynamical instability known as the `wiggle instability' and the thermal instability. The observed asymmetry must be transient: observations made tens of megayears in the past or in the future would often show an asymmetry in the opposite sense. Fluctuations of amplitude comparable to the observed asymmetry occur for a large fraction of time in our simulations, and suggest that the present is not an exceptional moment in the life of our Galaxy.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 475, Issue 2, p.2383-2402 (2018). A copy of the preprint is also available on the astro-ph server (arXiv:1707.03650).

The 30 Doradus nebula in the Large Magellanic Cloud (LMC) contains the massive starburst cluster NGC 2070 with a massive and probably younger stellar sub clump at its centre: R136. It is not clear how such a massive inner cluster could form several million years after the older stars in NGC 2070, given that stellar feedback is usually thought to expel gas and inhibit further star formation. Using the recently developed 1D feedback scheme WARPFIELD to scan a large range of cloud and cluster properties, we show that an age offset of several million years between the stellar populations is in fact to be expected given the interplay between feedback and gravity in a giant molecular cloud with a density $\geq$ 500 cm$^{-3}$ due to re-accretion of gas on to the older stellar population. Neither capture of field stars nor gas retention inside the cluster have to be invoked in order to explain the observed age offset in NGC 2070 as well as the structure of the interstellar medium around it.

This paper has been published in the Monthly Notices of the Royal Astronomical Society: Letters, Volume 473, Issue 1, p.L11-L15 (2018). A copy of the preprint is also available on the astro-ph server (arXiv:1710.02747).

Star clusters interact with the interstellar medium (ISM) in various ways, most importantly in the destruction of molecular star-forming clouds, resulting in inefficient star formation on galactic scales. On cloud scales, ionizing radiation creates H II regions, while stellar winds and supernovae (SNe) drive the ISM into thin shells. These shells are accelerated by the combined effect of winds, radiation pressure, and SN explosions, and slowed down by gravity. Since radiative and mechanical feedback is highly interconnected, they must be taken into account in a self-consistent and combined manner, including the coupling of radiation and matter. We present a new semi-analytic 1D feedback model for isolated massive clouds ($\geq$105 M$_{\odot}$) to calculate shell dynamics and shell structure simultaneously. It allows us to scan a large range of physical parameters (gas density, star formation efficiency, and metallicity) and to estimate escape fractions of ionizing radiation f$_{esc, I}$, the minimum star formation efficiency $\epsilon_{min}$ required to drive an outflow, and recollapse time-scales for clouds that are not destroyed by feedback. Our results show that there is no simple answer to the question of what dominates cloud dynamics, and that each feedback process significantly influences the efficiency of the others. We find that variations in natal cloud density can very easily explain differences between dense-bound and diffuse-open star clusters. We also predict, as a consequence of feedback, a 4-6 Myr age difference for massive clusters with multiple generations.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 470, Issue 4, p.4453-4472 (2017). A copy of the preprint is also available on the astro-ph server (arXiv:1704.04240).

Hydrodynamical simulations of star formation often do not possess the dynamic range needed to fully resolve the build-up of individual stars and star clusters, and thus have to resort to sub-grid models. A popular way to do this is by introducing Lagrangian sink particles, which replace contracting high-density regions at the point where the resolution limit is reached. A common problem then is how to assign fundamental stellar properties to sink particles, such as the distribution of stellar masses. We present a new and simple statistical method to assign stellar contents to sink particles. Once the stellar content is specified, it can be used to determine a sink particle's radiative output, supernovae rate or other feedback parameters that may be required in the calculations. Advantages of our method are: (I) it is simple to implement; (II) it guarantees that the obtained stellar populations are good samples of the initial mass function; (III) it can easily deal with infalling mass accreted at later times; and (IV) it does not put restrictions on the sink particles' masses in order to be used. The method works very well for sink particles that represent large star clusters and for which the stellar mass function is well sampled, but can also handle the transition to sink particles that represent a small number of stars.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 466, Issue 1, p.407-412 (2017). A copy of the preprint is also available on the astro-ph server (arXiv:1610.02538).

Aims: Polarization measurements of dust grains aligned with the magnetic field direction are an established technique for tracing large-scale field structures. In this paper we present a case study to investigate the conditions that need to be met to detect a characteristic magnetic field substructure that is embedded in such a large-scale field. A helical magnetic field with a surrounding hourglass-shaped field is expected from theoretical predictions and self-consistent magnetohydrodynamical (MHD) simulations to be present in the specific case of protostellar outflows. Hence, such an outflow environment is the perfect environment for our study.
Methods: We present synthetic polarization maps in the infrared and millimeter regime of simulations of protostellar outflows. The simulations were performed with the newly developed radiative transfer and polarization code POLARIS. The code is the first to include a self-consistent description of various alignment mechanisms such as the imperfect Davis-Greenstein (IDG) and the radiative torque (RAT) alignment. We investigated the effects of the grain size distribution, inclination, and applied alignment mechanism.
Results: We find that the IDG mechanism cannot produce any measurable polarization degree ($\geq$1%), whereas the RAT alignment produced polarization degrees of a few percent. Furthermore, we developed a method for identifying the origin of the polarization. We show that the helical magnetic field in the outflow can only be observed close to the outflow axis and at its tip, whereas in the surrounding regions the hourglass field in the foreground dominates the polarization. Furthermore, the polarization degree in the outflow lobe is lower than in the surroundings, in agreement with observations. We also find that the orientation of the polarization vector flips around at about a few hundred micrometers because of the transition from dichroic extinction to thermal re-emission. In order to avoid ambiguities when interpreting polarization data, we therefore suggest to observe in the far-infrared and millimeter regime. The actual grain size distribution has only little effect on the emerging polarization maps. Finally, we show that it is possible to observe the polarized radiation emerging from protostellar outflows with ALMA.

This paper has been published in the Astronomy & Astrophysics, Volume 603, id.A71, 14 pp.. A copy of the preprint is also available on the astro-ph server (arXiv:1703.02932).

We present three-dimensional magneto-hydrodynamical simulations of the self-gravitating interstellar medium (ISM) in a periodic (256 pc)$^3$ box with a mean number density of 0.5 cm$^{-3}$. At a fixed supernova rate we investigate the multi-phase ISM structure, H$_{2}$ molecule formation and density-magnetic field scaling for varying initial magnetic field strengths (0, $6\times 10^{-3}$, 0.3, 3 $\mu$G). All magnetic runs saturate at mass weighted field strengths of $\sim$ 1 $-$ 3 $\mu$G but the ISM structure is notably different. With increasing initial field strengths (from $6\times 10^{-3}$ to 3 $\mu$G) the simulations develop an ISM with a more homogeneous density and temperature structure, with increasing mass (from 5% to 85%) and volume filling fractions (from 4% to 85%) of warm (300 K $<$ T $<$ 8000 K) gas, with decreasing volume filling fractions (VFF) from $\sim$ 35% to $\sim$ 12% of hot gas (T $> 10^5$ K) and with a decreasing H$_{2}$ mass fraction (from 70% to $<$ 1%). Meanwhile the mass fraction of gas in which the magnetic pressure dominates over the thermal pressure increases by a factor of 10, from 0.07 for an initial field of $6\times 10^{-3}$ $\mu$G to 0.7 for a 3 $\mu$G initial field. In all but the simulations with the highest initial field strength self-gravity promotes the formation of dense gas and H$_{2}$, but does not change any other trends. We conclude that magnetic fields have a significant impact on the multi-phase, chemical and thermal structure of the ISM and discuss potential implications and limitations of the model.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 465, Issue 4, p.4611-4633. A copy of the preprint is also available on the astro-ph server (arXiv:1611.00585).

Supernova (SN) blast waves inject energy and momentum into the interstellar medium (ISM), control its turbulent multiphase structure and the launching of galactic outflows. Accurate modelling of the blast wave evolution is therefore essential for ISM and galaxy formation simulations. We present an efficient method to compute the input of momentum, thermal energy, and the velocity distribution of the shock-accelerated gas for ambient media (densities of 0.1 $\geq$ n$_0$ [cm$^{-3}$] $\geq$ 100) with uniform (and with stellar wind blown bubbles), power-law, and turbulent (Mach numbers M from 1 - 100) density distributions. Assuming solar metallicity cooling, the blast wave evolution is followed to the beginning of the momentum conserving snowplough phase. The model recovers previous results for uniform ambient media. The momentum injection in wind-blown bubbles depend on the swept-up mass and the efficiency of cooling, when the blast wave hits the wind shell. For power-law density distributions with n(r) $\sim$ r$^{-2}$ (for n(r) $\gt$ n$_{floor}$) the amount of momentum injection is solely regulated by the background density n$_{floor}$ and compares to $n_{uni} = n_{floor}$. However, in turbulent ambient media with lognormal density distributions the momentum input can increase by a factor of 2 (compared to the homogeneous case) for high Mach numbers. The average momentum boost can be approximated as $p_{turb}/p_{0} =23.07 (\frac{n_{0,turb}}{1 cm^{-3}})^{-0.12} + 0.82 (ln (1+b^2M^2))^{1.49}(\frac{n_{0,turb}}{1 cm^{-3}})^{-1.6}$. The velocity distributions are broad as gas can be accelerated to high velocities in low-density channels. The model values agree with results from recent, computationally expensive, three-dimensional simulations of SN explosions in turbulent media.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 460, Issue 3, p.2962-2978 (2016). A copy of the preprint is also available on the astro-ph server (arXiv:1604.04395).

Two-thirds of long duration gamma-ray bursts (GRBs) show soft X-ray absorption in excess of the Milky Way. The column densities of metals inferred from UV and optical spectra differ from those derived from soft X-ray spectra, at times by an order of magnitude, with the latter being higher. The origin of the soft X-ray absorption excess observed in GRB X-ray afterglow spectra remains a heavily debated issue, which has resulted in numerous investigations on the effect of hot material both internal and external to the GRB host galaxy on our X-ray afterglow observations. Nevertheless, all models proposed so far have either only been able to account for a subset of our observations (i.e. at $z > 2$), or they have required fairly extreme conditions to be present within the absorbing material. In this paper, we investigate the absorption of the GRB afterglow by a collisionally ionised and turbulent interstellar medium (ISM). We find that a dense (3 per cubic centimeters) collisionally ionised ISM could produce UV/optical and soft X-ray absorbing column densities that differ by a factor of 10, however the UV/optical and soft X-ray absorbing column densities for such sightlines and are 2-3 orders of magnitude lower in comparison to the GRB afterglow spectra. For those GRBs with a larger soft X-ray excess of up to an order of magnitude, the contribution in absorption from a turbulent ISM as considered here would ease the required conditions of additional absorbing components, such as the GRB circumburst medium and intergalactic medium.

This paper has been published in the Astronomy & Astrophysics, Volume 595, id.A24, 10. A copy of the preprint is also available on the astro-ph server (arXiv:1603.05123).

The halo of the Milky Way contains a hot plasma with a surface brightness in soft X-rays of the order $10^{-12}{\rm~erg} {\rm~cm}^{-2} {\rm~s}^{-1} {\rm~deg}^{-2}$. The origin of this gas is unclear, but so far numerical models of galactic star formation have failed to reproduce such a large surface brightness by several orders of magnitude. In this paper, we analyze simulations of the turbulent, magnetized, multi-phase interstellar medium including thermal feedback by supernova explosions as well as cosmic-ray feedback. We include a time-dependent chemical network, self-shielding by gas and dust, and self-gravity. Pure thermal feedback alone is sufficient to produce the observed surface brightness, although it is very sensitive to the supernova rate. Cosmic rays suppress this sensitivity and reduce the surface brightness because they drive cooler outflows. Self-gravity has by far the largest effect because it accumulates the diffuse gas in the disk in dense clumps and filaments, so that supernovae exploding in voids can eject a large amount of hot gas into the halo. This can boost the surface brightness by several orders of magnitude. Although our simulations do not reach a steady state, all simulations produce surface brightness values of the same order of magnitude as the observations, with the exact value depending sensitively on the simulation parameters. We conclude that star formation feedback alone is sufficient to explain the origin of the hot halo gas, but measurements of the surface brightness alone do not provide useful diagnostics for the study of galactic star formation.

This paper has been published in the The Astrophysical Journal Letters, Volume 813, Issue 2, article id. L27, 7 pp. (2015). A copy of the preprint is also available on the astro-ph server (arXiv:1510.06563).

We present three-dimensional magneto-hydrodynamical simulations of the self-gravitating interstellar medium (ISM) in a periodic (256 pc)$^3$ box with a mean number density of 0.5 cm$^{-3}$. At a fixed supernova rate we investigate the multi-phase ISM structure, H$_{2}$ molecule formation and density-magnetic field scaling for varying initial magnetic field strengths (0, $6\times 10^{-3}$, 0.3, 3 $\mu$G). All magnetic runs saturate at mass weighted field strengths of $\sim$ 1 $-$ 3 $\mu$G but the ISM structure is notably different. With increasing initial field strengths (from $6\times 10^{-3}$ to 3 $\mu$G) the simulations develop an ISM with a more homogeneous density and temperature structure, with increasing mass (from 5% to 85%) and volume filling fractions (from 4% to 85%) of warm (300 K $<$ T $<$ 8000 K) gas, with decreasing volume filling fractions (VFF) from $\sim$ 35% to $\sim$ 12% of hot gas (T $> 10^5$ K) and with a decreasing H$_{2}$ mass fraction (from 70% to $<$ 1%). Meanwhile the mass fraction of gas in which the magnetic pressure dominates over the thermal pressure increases by a factor of 10, from 0.07 for an initial field of $6\times 10^{-3}$ $\mu$G to 0.7 for a 3 $\mu$G initial field. In all but the simulations with the highest initial field strength self-gravity promotes the formation of dense gas and H$_{2}$, but does not change any other trends. We conclude that magnetic fields have a significant impact on the multi-phase, chemical and thermal structure of the ISM and discuss potential implications and limitations of the model.

This paper has been published in the Monthly Notices of the Royal Astronomical Society, Volume 449, Issue 1, p.1057-1075. A copy of the preprint is also available on the astro-ph server (arXiv:1411.0009).