(Context) Monte Carlo radiative transfer (MCRT) is a widely used technique to model the interaction between radiation and a medium, and plays an important role in astrophysical modelling and when comparing those models with observations. (Aims) In this work, we present a novel approach to MCRT that addresses the challenging memory access patterns of traditional MCRT algorithms, which hinder optimal performance of MCRT simulations on modern hardware with a complex memory architecture. (Methods) We reformulate the MCRT photon packet life cycle as a task-based algorithm, whereby the computation is broken down into small tasks that are executed concurrently. Photon packets are stored in intermediate buffers, and tasks propagate photon packets through small parts of the computational domain, moving them from one buffer to another in the process. (Results) Using the implementation of the new algorithm in the photoionization MCRT code CMacIonize 2.0, we show that the decomposition of the MCRT grid into small parts leads to a significant performance gain during the photon packet propagation phase, which constitutes the bulk of an MCRT algorithm, as a result of better usage of memory caches. Our new algorithm is a factor 2 to 4 faster than an equivalent traditional algorithm and shows good strong scaling up to 30 threads. We briefly discuss how our new algorithm could be adjusted or extended to other astrophysical MCRT applications. (Conclusions) We show that optimising the memory access patterns of a memory-bound algorithm such as MCRT can yield significant performance gains.

We present the new open source C++-based Python library CosTuuM that can be used to generate infrared absorption and emission coefficients for arbitrary mixtures of spheroidal dust grains that are (partially) aligned with a magnetic field. We outline the algorithms underlying the software, demonstrate the accuracy of our results using benchmarks from literature, and use our tool to investigate some commonly used approximative recipes. We find that the linear polarization fraction for a partially aligned dust grain mixture can be accurately represented by an appropriate linear combination of perfectly aligned grains and grains that are randomly oriented, but that the commonly used picket fence alignment breaks down for short wavelengths. We also find that for a fixed dust grain size, the absorption coefficients and linear polarization fraction for a realistic mixture of grains with various shapes cannot both be accurately represented by a single representative grain with a fixed shape, but that instead an average over an appropriate shape distribution should be used. Insufficient knowledge of an appropriate shape distribution is the main obstacle in obtaining accurate optical properties. CosTuuM is available as a standalone Python library and can be used to generate optical properties to be used in radiative transfer applications.

We present infrared luminosity functions and dust mass functions for the EAGLE cosmological simulation, based on synthetic multi-wavelength observations generated with the SKIRT radiative transfer code. In the local Universe, we reproduce the observed infrared luminosity and dust mass functions very well. Some minor discrepancies are encountered, mainly in the high luminosity regime, where the EAGLE-SKIRT luminosity functions mildly but systematically underestimate the observed ones. The agreement between the EAGLE-SKIRT infrared luminosity functions and the observed ones gradually worsens with increasing lookback time. Fitting modified Schechter functions to the EAGLE-SKIRT luminosity and dust mass functions at different redshifts up to $z=1$, we find that the evolution is compatible with pure luminosity/mass evolution. The evolution is relatively mild: within this redshift range, we find an evolution of $L_{\star,250}\propto(1+z)^{1.68}$, $L_{\star,\text{TIR}}\propto(1+z)^{2.51}$ and $M_{\star,\text{dust}}\propto(1+z)^{0.83}$ for the characteristic luminosity/mass. For the luminosity/mass density we find $\varepsilon_{250}\propto(1+z)^{1.62}$, $\varepsilon_{\text{TIR}}\propto(1+z)^{2.35}$ and $\rho_{\text{dust}}\propto(1+z)^{0.80}$, respectively. The mild evolution of the dust mass density is in relatively good agreement with observations, but the slow evolution of the infrared luminosity underestimates the observed luminosity evolution significantly. We argue that these differences can be attributed to increasing limitations in the radiative transfer treatment due to increasingly poorer resolution, combined with a slower than observed evolution of the SFR density in the EAGLE simulation and the lack of AGN emission in our EAGLE-SKIRT post-processing recipe.

There is strong evidence that the diffuse ionized gas (DIG) in disc galaxies is photoionized by radiation from UV luminous O and B stars in the galactic disc, both from observations and detailed numerical models. However, it is still not clear what mechanism is responsible for providing the necessary pressure support for a diffuse gas layer at kpc-scale above the disc. In this work we investigate if the pressure increase caused by photoionization can provide this support. We run self-consistent radiation hydrodynamics models of a gaseous disc in an external potential. We find that photoionization feedback can drive low levels of turbulence in the dense galactic disc, and that it provides pressure support for an extended diffuse gas layer. Our results show that there is a natural fine-tuning between the total ionizing radiation budget of the sources in the galaxy and the amount of gas in the different ionization phases of the ISM, and provide the first fully consistent radiation hydrodynamics model of the DIG.

The formation of massive stars is a long standing problem. Although a number of theories of massive star formation exist, ideas appear to converge to a disk-mediated accretion scenario. Here we present radiative hydrodynamic simulations of a star accreting mass via a disk embedded in a torus. We use a Monte Carlo based radiation hydrodynamics code to investigate the impact that ionizing radiation has on the torus. Ionized regions in the torus midplane are found to be either gravitationally trapped or in pressure driven expansion depending on whether or not the size of the ionized region exceeds a critical radius. Trapped Hii regions in the torus plane allow accretion to progress, while expanding Hii regions disrupt the accretion torus preventing the central star from aggregating more mass, thereby setting the star's final mass. We obtain constraints for the luminosities and torus densities that lead to both scenarios.

This paper investigates the gravitational trapping of HII regions predicted by steady-state analysis using radiation hydrodynamical simulations. We present idealised spherically symmetric radiation hydrodynamical simulations of the early evolution of HII regions including the gravity of the central source. As with analytic steady state solutions of spherically symmetric ionised Bondi accretion flows, we find gravitationally trapped HII regions with accretion through the ionisation front onto the source. We found that, for a constant ionising luminosity, fluctuations in the ionisation front are unstable. This instability only occurs in this spherically symmetric accretion geometry. In the context of massive star formation, the ionising luminosity increases with time as the source accretes mass. The maximum radius of the recurring HII region increases on the accretion timescale until it reaches the sonic radius, where the infall velocity equals the sound speed of the ionised gas, after which it enters a pressure-driven expansion phase. This expansion prevents accretion of gas through the ionisation front, the accretion rate onto the star decreases to zero, and it stops growing from accretion. Because of the time required for any significant change in stellar mass and luminosity through accretion our simulations keep both mass and luminosity constant and follow the evolution from trapped to expanding in a piecewise manner. Implications of this evolution of HII regions include a continuation of accretion of material onto forming stars for a period after the star starts to emit ionising radiation, and an extension of the lifetime of ultracompact HII regions.

We investigate the general stability of 1D spherically symmetric ionized Bondi accretion onto a massive object in the specific context of accretion onto a young stellar object. We first derive a new analytic expression for a steady state two temperature solution that predicts the existence of compact and hypercompact HII regions. We then show that this solution is only marginally stable if ionization is treated self-consistently. This leads to a recurring collapse of the HII region over time. We derive a semi-analytic model to explain this instability, and test it using spatially converged 1D radiation hydrodynamical simulations. We discuss the implications of the 1D instability on 3D radiation hydrodynamics simulations of supersonic accreting flows.

We present the public Monte Carlo photoionization and moving-mesh radiation hydrodynamics code CMacIonize, which can be used to simulate the self-consistent evolution of HII regions surrounding young O and B stars, or other sources of ionizing radiation. The code combines a Monte Carlo photoionization algorithm that uses a complex mix of hydrogen, helium and several coolants in order to self-consistently solve for the ionization and temperature balance at any given type, with a standard first order hydrodynamics scheme. The code can be run as a post-processing tool to get the line emission from an existing simulation snapshot, but can also be used to run full radiation hydrodynamical simulations. Both the radiation transfer and the hydrodynamics are implemented in a general way that is independent of the grid structure that is used to discretize the system, allowing it to be run both as a standard fixed grid code, but also as a moving-mesh code.

The large vertical scale heights of the diffuse ionised gas (DIG) in disc galaxies are challenging to model, as hydrodynamical models including only thermal feedback seem to be unable to support gas at these heights. In this paper, we use a three dimensional Monte Carlo radiation transfer code to post-process disc simulations of the Simulating the Life-Cycle of Molecular Clouds (SILCC) project that include feedback by cosmic rays. We show that the more extended discs in simulations including cosmic ray feedback naturally lead to larger scale heights for the DIG which are more in line with observed scale heights. We also show that including a fiducial cosmic ray heating term in our model can help to increase the temperature as a function of disc scale height, but fails to reproduce observed DIG nitrogen and sulphur forbidden line intensities. We show that, to reproduce these line emissions, we require a heating mechanism that affects gas over a larger density range than is achieved by cosmic ray heating, which can be achieved by fine tuning the total luminosity of ionising sources to get an appropriate ionising spectrum as a function of scale height. This result sheds a new light on the relation between forbidden line emissions and temperature profiles for realistic DIG gas distributions.

We introduce the moving mesh code Shadowfax, which can be used to evolve a mixture of gas, subject to the laws of hydrodynamics and gravity, and any collisionless fluid only subject to gravity, such as cold dark matter or stars. The code is written in C++ and its source code is made available to the scientific community under the GNU Affero General Public License. We outline the algorithm and the design of our implementation, and demonstrate its validity through the results of a set of basic test problems, which are also part of the public version. We also compare Shadowfax with a number of other publicly available codes using different hydrodynamical integration schemes, illustrating the advantages and disadvantages of the moving mesh technique.

Simulating dwarf galaxy halos in a reionizing Universe puts severe constraints on the sub-grid model employed in the simulations. Using the same sub-grid model that works for simulations without a UV-background (UVB) results in gas poor galaxies that stop forming stars very early on, except for halos with high masses. This is in strong disagreement with observed galaxies, which are gas rich and star forming down to a much lower mass range. To resolve this discrepancy, we ran a large suite of isolated dwarf galaxy simulations to explore a wide variety of sub-grid models and parameters, including timing and strength of the UVB, strength of the stellar feedback, and metallicity dependent Pop III feedback. We compared these simulations to observed dwarf galaxies by means of the baryonic Tully-Fisher relation (BTFR), which links the baryonic content of a galaxy to the observationally determined strength of its gravitational potential. We found that the results are robust to changes in the UVB. The strength of the stellar feedback shifts the results on the BTFR, but does not help to form gas rich galaxies at late redshifts. Only by including Pop III feedback are we able to produce galaxies that lie on the observational BTFR and that have neutral gas and ongoing star formation at redshift zero.

Low-mass dwarf galaxies are very sensitive test-beds for theories of cosmic structure formation since their weak gravitational fields allow the effects of the relevant physical processes to clearly stand out. Up to now, no unified account exists of the sometimes seemingly conflicting properties of the faintest isolated dwarfs in and around the Local Group, such as Leo T and the recently discovered Leo P and Pisces A systems. Using new numerical simulations, we show that this serious challenge to our understanding of galaxy formation can be effectively resolved by taking into account the regulating influence of the ultraviolet radiation of the first population of stars on a dwarf's star formation rate while otherwise staying within the standard cosmological paradigm for structure formation. These simulations produce faint, gas-dominated, star-forming dwarf galaxies that lie on the baryonic Tully-Fisher relation and that successfully reproduce a broad range of chemical, kinematical, and structural observables of real late-type dwarf galaxies. Furthermore, we stress the importance of obtaining properties of simulated galaxies in a manner as close as possible to the typically employed observational techniques.

We investigate the evolution of dwarf galaxies using N -body/SPH simulations that incorporate their formation histories through merger trees constructed using the ex- tended Press-Schechter formalism. The simulations are computationally cheap and have high spatial resolution. We compare the properties of galaxies with equal final mass but with different merger histories with each other and with those of observed dwarf spheroidals and irregulars. We show that the merger history influences many observable dwarf galaxy proper- ties. We identify two extreme cases that make this influence stand out most clearly: (i) merger trees with one massive progenitor that grows through relatively few mergers and (ii) merger trees with many small progenitors that merge only quite late. At a fixed halo mass, a type (i) tree tends to produce galaxies with larger stellar masses, larger half-light radii, lower central surface brightness, and since fewer potentially an- gular momentum cancelling mergers are required to build up the final galaxy, a higher specific angular momentum, compared with a type (ii) tree. We do not perform full-fledged cosmological simulations and therefore cannot hope to reproduce all observed properties of dwarf galaxies. However, we show that the simulated dwarfs are not unsimilar to real ones.

Using computer simulations, we explored gaseous infall as a possible explanation for the starburst phase in Blue Compact Dwarf galaxies. We simulate a cloud impact by merging a spherical gas cloud into an isolated dwarf galaxy. We investigated which conditions were favourable for triggering a burst and found that the orbit and the mass of the gas cloud play an important role. We discuss the metallicity, the kinematical properties, the internal dynamics and the gas, stellar and dark matter distribution of the simulations during a starburst. We find that these are in good agreement with observations and depending on the set-up (e.g. rotation of the host galaxy, radius of the gas cloud), our bursting galaxies can have qualitatively very different properties. Our simulations offer insight in how starbursts start and evolve. Based on this, we propose what postburst dwarf galaxies will look like.

We present a detailed analysis of the formation, evolution, and possible longevity of metallicity gradients in simulated dwarf galaxies. Specifically, we investigate the role of potentially orbit-changing processes such as radial stellar migration and dynamical heating in shaping or destroying these gradients. We also consider the influence of the star formation density threshold, investigating both a low and high value (0.1 amu/$cm^{3}$ - 100 amu/$cm^{3}$). The Nbody-SPH models that we use to self-consistently form and evolve dwarf galaxies in isolation show that, in the absence of significant angular momentum, metallicity gradients are gradually built up during the evolution of the dwarf galaxy, by ever more centrally concentrated star formation adding to the overall gradient. Once formed, they are robust, survive easily in the absence of external disturbances and hardly decline over several Gyr, and they agree well with observed metallicity gradients of dwarf galaxies in the Local Group. The underlying orbital displacement of stars is quite limited in our models, being of the order of only fractions of the half light radius over time-spans of 5 to 10 Gyr in all star formation schemes. This is contrary to the strong radial migration found in massive disc galaxies, caused by scattering of stars off the corotation resonance of spiral structures. In the dwarf regime the stellar body only seems to undergo mild dynamical heating, due to the lack of long-lived spiral structures and/or discs. The density threshold has profound influences on the star formation mode of the models, but only a minor influence on the evolution of metallicity gradients. Increasing the threshold 1000-fold causes comparatively stronger dynamical heating of the stellar body due to the increased turbulence and the scattering of stars off dense gas clouds, but the effect remains very limited in absolute terms.

In this paper, we present a new calculation of composition-dependent radiative cooling and heating curves of low-density gas, intended primarily for use in numerical simulations of galaxy formation and evolution. These curves depend on only five parameters: temperature, density, redshift, [Fe/H], and [Mg/Fe]. They are easily tabulated and can be efficiently interpolated during a simulation. The ionization equilibrium of 14 key elements is determined for temperatures between 10K and 10^9K and densities up to 100 amu/cm^3 taking into account collisional and radiative ionization, by the cosmic UV background and an interstellar radiation field, and by charge-transfer reactions. These elements, ranging from H to Ni, are the ones most abundantly produced and/or released by SNIa, SNII, and intermediate-mass stars. Self-shielding of the gas at high densities by neutral Hydrogen is taken into account in an approximate way by exponentially suppressing the H-ionizing part of the cosmic UV background for HI densities above a threshold density of n_HI,crit=0.007 cm^-3. We discuss how the ionization equilibrium, and the cooling and heating curves depend on the physical properties of the gas. The main advantage of the work presented here is that, within the confines of a well-defined chemical evolution model and adopting the ionization equilibrium approximation, it provides accurate cooling and heating curves for a wide range of physical and chemical gas properties, including the effects of self-shielding. The latter is key to resolving the formation of cold, neutral, high-density clouds suitable for star formation in galaxy simulations.

Simulation codes for galaxy formation and evolution take on board as many physical processes as possible beyond the standard gravitational and hydrodynamical physics. Most of this extra physics takes place below the resolution level of the simulations and is added in a sub-grid fashion. However, these sub-grid processes affect the macroscopic hydrodynamical properties of the gas and thus couple to the on-grid physics that is explicitly integrated during the simulation. In this paper, we focus on the link between partial ionization and the hydrodynamical equations. We show that the energy stored in ions and free electrons constitutes a potential energy term which breaks the linear dependence of the internal energy on temperature. Correctly taking into account ionization hence requires modifying both the equation of state and the energy-temperature relation. We implemented these changes in the cosmological simulation code Gadget2. As an example of the effects of these changes, we study the propagation of Sedov-Taylor shock waves through an ionizing medium. This serves as a proxy for the absorption of supernova feedback energy by the interstellar medium. Depending on the density and temperature of the surrounding gas, we find that up to 50% of the feedback energy is spent ionizing the gas rather than heating it. Thus, it can be expected that properly taking into account ionization effects in galaxy evolution simulations will drastically reduce the effects of thermal feedback. To the best of our knowledge, this potential energy term is not used in current simulations of galaxy formation and evolution.