Search

Using the First Light And Reionisation Epoch Simulations (\textsc{Flares}) we explore the dust driven properties of massive high-redshift galaxies at $z\in[5,10]$. By post-processing the galaxy sample using the radiative transfer code \textsc{skirt} we obtain the full spectral energy distribution. We explore the resultant luminosity functions, IRX-$\beta$ relations as well as the luminosity-weighted dust temperatures in the Epoch of Reionisation (EoR). We find that most of our results are in agreement with the current set of observations, but under-predict the number densities of bright IR galaxies, which are extremely biased towards the most overdense regions. We see that the \textsc{Flares} IRX-$\beta$ relation (for $5\le z\le8$) predominantly follows the local starburst relation. The IRX shows an increase with stellar mass, plateauing at the high-mass end ($\sim10^{10}$M$_{\odot}$) and shows no evolution in the median normalisation with redshift. We also look at the dependence of the peak dust temperature ($T_{\mathrm{peak}}$) on various galaxy properties including the stellar mass, IR luminosity and sSFR, finding the correlation to be strongest with sSFR. The luminosity-weighted dust temperatures increase towards higher redshifts, with the slope of the $T_{\mathrm{peak}}$ - redshift relation showing a higher slope than the lower redshift relations obtained from previous observational and theoretical works. The results from \textsc{Flares}, which is able to provide a better statistical sample of high-redshift galaxies compared to other simulations, provides a distinct vantage point for the high-redshift Universe., Comment: 19 pages, 16 figures, Published in MNRAS. Data link: https://flaresimulations.github.io/data.html

Dust emission from high-redshift galaxies gives us a clue to the origin and evolution of dust in the early Universe. Previous studies have shown that different sources of dust (stellar dust production and dust growth in dense clouds) predict different ultraviolet (UV) extinction curves for galaxies at $z\sim 7$ but that the observed attenuation curves depend strongly on the geometry of dust and star distributions. Thus, we perform radiative transfer calculations under the dust-stars geometries computed by a cosmological hydrodynamic simulation (IllustrisTNG). This serves to investigate the dust attenuation curves predicted from `realistic' geometries. We choose objects with stellar mass and star formation rate appropriate for Lyman break galaxies at $z\sim 7$. We find that the attenuation curves are very different from the original extinction curves in most of the galaxies. This makes it difficult to constrain the dominant dust sources from the observed attenuation curves. We further include infrared dust emission in the analysis and plot the infrared excess (IRX)-UV spectral slope ($\beta$) diagram. We find that different sources of dust cause different IRX-$\beta$ relations for the simulated galaxies. In particular, if dust growth is the main source of dust, a variation of dust-to-metal ratio causes a more extended sequence with smaller IRX in the IRX-$\beta$ diagram. Thus, the comprehensive analysis of the abundances of dust and metals, the UV slope, and the dust emission could provide a clue to the dominant dust sources in the Universe., Comment: 11 pages, 6 figures, accepted for publication in MNRAS

Cold gas in galaxies provides a crucial test to evaluate the realism of cosmological hydrodynamical simulations. To extract the atomic and molecular hydrogen properties of the simulated galaxy population, postprocessing methods taking the local UV field into account are required. We improve upon previous studies by calculating realistic UV fields with the dust radiative transfer code SKIRT to model the atomic-to-molecular transition in TNG50, the highest-resolution run of the IllustrisTNG suite. Comparing integrated quantities such as the HI mass function, we study to what detail the UV field needs to be modelled in order to calculate realistic cold gas properties. We then evaluate new, spatially resolved comparisons for cold gas in galaxies by exploring synthetic maps of atomic hydrogen at redshift zero and compare them to 21-cm observations of local galaxies from the WHISP survey. In terms of non-parametric morphologies, we find that TNG50 HI maps are less concentrated than their WHISP counterparts (median $\Delta C\approx0.3$), due in part to central HI deficits related to the ejective character of supermassive black hole feedback in TNG. In terms of the HI column density distribution function, we find discrepancies between WHISP and IllustrisTNG that depend on the total HI abundance in these datasets as well as the postprocessing method. To fully exploit the synergy between cosmological simulations and upcoming deep HI/H2 data, we advocate the use of accurate methods to estimate the UV radiation field and to generate mock maps., Comment: 21 pages, 10 figures (main text). Accepted to MNRAS

Dust has been detected in high-redshift ($z>5$) galaxies but its origin is still being debated. Dust production in high-redshift galaxies could be dominated by stellar production or by accretion (dust growth) in the interstellar medium. Previous studies have shown that these two dust sources predict different grain size distributions, which lead to significantly different extinction curves. In this paper, we investigate how the difference in the extinction curves affects the dust attenuation properties of galaxies by performing radiative transfer calculations. To examine the major effects of the dust--stars distribution geometry, we adopt two representative cases in spherical symmetry: the well-mixed geometry (stars and dust are homogeneously mixed) and the two-layer geometry (young stars are more concentrated in the centre). In both cases, we confirm that the attenuation curve can be drastically steepened by scattering and by different optical depths between young and old stellar populations, and can be flattened by the existence of unobscured stellar populations. We can reproduce similar attenuation curves even with very different extinction curves. Thus, we conclude that it is difficult to distinguish the dust sources only with attenuation curves. However, if we include information on dust emission and plot the IRX (infrared excess)--$\beta$ (ultraviolet spectral slope) relation, different dust sources predict different positions in the IRX--$\beta$ diagram. A larger $\beta$ is preferred under a similar IRX if dust growth is the dominant dust source., Comment: 11 pages, 4 figures, accepted for publication in MNRAS

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., Comment: 13 pages, 6 pages, accepted for publication in MNRAS

Context: The Monte Carlo method is probably the most widely used approach to solve the radiative transfer problem, especially in a general 3D geometry. The physical processes of emission, absorption, and scattering are easily incorporated in the Monte Carlo framework. Net stimulated emission, or absorption with a negative cross section, does not fit this method, however. Aims: We explore alterations to the standard photon packet life cycle in Monte Carlo radiative transfer that allow the treatment of net stimulated emission without loss of generality or efficiency. Methods: We present the explicit absorption technique that allows net stimulated emission to be handled efficiently. It uses the scattering rather than the extinction optical depth along a photon packet's path to randomly select the next interaction location, and offers a separate, deterministic treatment of absorption. We implemented the technique in a special-purpose Monte Carlo code for a two-stream 1D radiative transfer problem and in the fully featured 3D code SKIRT, and we studied its overall performance using quantitative statistical tests. Results: Our special-purpose code is capable of recovering the analytical solutions to the two-stream problem in all regimes, including the one of strong net stimulated emission. The implementation in SKIRT is straightforward, as the explicit absorption technique easily combines with the variance reduction and acceleration techniques already incorporated. In general, explicit absorption tends to improve the efficiency of the Monte Carlo routine in the regime of net absorption. Conclusions: Explicit absorption allows the treatment of net stimulated emission in Monte Carlo radiative transfer, it interfaces smoothly with... (abridged), Comment: 14 pages, accepted for publication in A&A

We apply the radiative transfer (RT) code SKIRT on a sample of ~14000 low-redshift (z, The full galaxy SEDs and the broadband fluxes at z

Context: Properly modelling scattering by interstellar dust grains requires a good characterisation of the scattering phase function. The Henyey-Greenstein phase function has become the standard for describing anisotropic scattering by dust grains, but it is a poor representation of the real scattering phase function outside the optical range. Aims: We investigate alternatives for the Henyey-Greenstein phase function that would allow the scattering properties of dust grains to be described. Our goal is to find a balance between realism and complexity: the scattering phase function should be flexible enough to provide an accurate fit to the scattering properties of dust grains over a wide wavelength range, and it should be simple enough to be easy to handle, especially in the context of radiative transfer calculations. Methods: We fit various analytical phase functions to the scattering phase function corresponding to the BARE-GR-S model, one of the most popular and commonly adopted models for interstellar dust. We weigh the accuracy of the fit against the number of free parameters in the analytical phase functions. Results: We confirm that the Henyey-Greenstein phase functions poorly describe scattering by dust grains, particularly at ultraviolet (UV) wavelengths, with relative differences of up to 50%. The Draine phase function alleviates this problem at near-infrared (NIR) wavelengths, but not in the UV. The two-term Reynolds-McCormick phase function, recently advocated in the context of light scattering in nanoscale materials and aquatic media, describes the BARE-GR-S data very well, but its five free parameters are degenerate. We propose a simpler phase function, the two-term ultraspherical-2 (TTU2) phase function, that also provides an excellent fit to the BARE-GR-S phase function over the entire UV-NIR wavelength range. (Abridged), Comment: 12 pages, 7 figures, accepted for publication in A&A

We present redshift-zero synthetic observational data considering dust attenuation and dust emission for the thirty galaxies of the Auriga project, calculated with the SKIRT radiative transfer code. The post-processing procedure includes components for star-forming regions, stellar sources, and diffuse dust taking into account stochastic heating of dust grains. This allows us to obtain realistic high-resolution broadband images and fluxes from ultraviolet to sub-millimeter wavelengths. For the diffuse dust component, we consider two mechanisms for assigning dust to gas cells in the simulation. In one case, only the densest or the coldest gas cells are allowed to have dust, while in the other case this condition is relaxed to allow a larger number of dust-containing cells. The latter approach yields galaxies with a larger radial dust extent and an enhanced dust presence in the inter-spiral regions. At a global scale, we compare Auriga galaxies with observations by deriving dust scaling relations using SED fitting. At a resolved scale, we make a multi-wavelength morphological comparison with nine well-resolved spiral galaxies from the DustPedia observational database. We find that for both dust assignment methods, although the Auriga galaxies show a good overall agreement with observational dust properties, they exhibit a slightly higher specific dust mass. The multi-wavelength morphological analysis reveals a good agreement between the Auriga and the observed galaxies in the optical wavelengths. In the mid and far-infrared wavelengths, Auriga galaxies appear smaller and more centrally concentrated in comparison to their observed counterparts. We publicly release the multi-observer images and fluxes in 50 commonly used broadband filters., 19 pages, 20 figures, accepted for publication in MNRAS

The cosmic spectral energy distribution (CSED) is the total emissivity as a function of wavelength of galaxies in a given cosmic volume. We compare the observed CSED from the UV to the submm to that computed from the EAGLE cosmological hydrodynamical simulation, post-processed with stellar population synthesis models and including dust radiative transfer using the SKIRT code. The agreement with the data is better than 0.15 dex over the entire wavelength range at redshift $z=0$, except at UV wavelengths where the EAGLE model overestimates the observed CSED by up to a factor 2. Global properties of the CSED as inferred from CIGALE fits, such as the stellar mass density, mean star formation density, and mean dust-to-stellar-mass ratio, agree to within better than 20 per cent. At higher redshift, EAGLE increasingly underestimates the CSED at optical-NIR wavelengths with the FIR/submm emissivity underestimated by more than a factor of 5 by redshift $z=1$. We believe that these differences are due to a combination of incompleteness of the EAGLE-SKIRT database, the small simulation volume and the consequent lack of luminous galaxies, and our lack of knowledge on the evolution of the characteristics of the interstellar dust in galaxies. The impressive agreement between the simulated and observed CSED at lower $z$ confirms that the combination of EAGLE and SKIRT dust processing yields a fairly realistic representation of the local Universe., 14 pages, 5 figures, accepted for publication in MNRAS

Galaxy kinematics and gravitational lensing are two complementary ways to constrain the distribution of dark matter on galaxy scales. The typical dark matter density profiles adopted in dynamical studies cannot easily be adopted in lensing studies. Ideally, a mass model should be used that has the global characteristics of realistic dark matter distributions, and that allows for an analytical calculation of the magnifications and deflection angles. A simple model with these properties, the broken-power-law (BPL) model, has very recently been introduced. We examine the dynamical structure of the family of BPL models. We derive simple closed expressions for basic dynamical properties, and study the distribution function under the assumption of velocity isotropy. We find that none of the BPL models with realistic parameters has an isotropic distribution function that is positive over the entire phase space, implying that the BPL models cannot be supported by an isotropic velocity distribution, or models with a more radially anisotropic orbital structure. This result limits the attractiveness of the BPL family as a tool for lensing studies to some degree. More generally, we find that not all members of the general family of double power-law or Zhao models, often used to model dark matter haloes, can be supported by an isotropic or radially anisotropic distribution function. In other words, the distribution function may become negative even for spherically symmetric models with a well-behaved density profile., 12 pages, 5 figures, accepted for publication in MNRAS

(Aims) In this work, we predict the far-infrared polarisation signal emitted by non-spherical dust grains in nearby galaxies. We determine the angular resolution and sensitivity required to study the magnetic field configuration in these galaxies. (Methods) We post-process a set of Milky Way like galaxies from the Auriga project, assuming a dust mix consisting of spheroidal dust grains that are partially aligned with the model magnetic field. We constrain our dust model using Planck 353 GHz observations of the Milky Way. This model is then extrapolated to shorter wavelengths that cover the peak of interstellar dust emission and to observations of arbitrarily oriented nearby Milky Way like galaxies. (Results) Assuming an intrinsic linear polarisation fraction that does not vary significantly with wavelength for wavelengths longer than 50 micron, we predict a linear polarisation fraction with a maximum of $10-15\%$ and a median value of $\approx{}7\%$ for face-on galaxies and $\approx{}3\%$ for edge-on galaxies. The polarisation fraction anti-correlates with the line of sight density and with the angular dispersion function which expresses the large scale order of the magnetic field perpendicular to the line of sight. The maximum linear polarisation fraction agrees well with the intrinsic properties of the dust model. The true magnetic field orientation can be traced along low density lines of sight when it is coherent along the line of sight. These results also hold for nearby galaxies, where a coherent magnetic field structure is recovered over a range of different broad bands. (Conclusions) Polarised emission from non-spherical dust grains accurately traces the large scale structure of the galactic magnetic field in Milky Way like galaxies, with expected maximum linear polarisation fractions of $10-15\%$. To resolve this maximum, a spatial resolution of at least 1 kpc is required., Comment: 22 pages, 19 figures, accepted for publication in Astronomy and Astrophysics

There is a severe tension between the observed star formation rate (SFR) - stellar mass (${\rm M}_{\star}$) relations reported by different authors at $z = 1-4$. In addition, the observations have not been successfully reproduced by state-of-the-art cosmological simulations which tend to predict a factor of 2-4 smaller SFRs at a fixed ${\rm M}_{\star}$. We examine the evolution of the SFR$-{\rm M}_{\star}$ relation of $z = 1-4 $ galaxies using the SKIRT simulated spectral energy distributions of galaxies sampled from the EAGLE simulations. We derive SFRs and stellar masses by mimicking different observational techniques. We find that the tension between observed and simulated SFR$-{\rm M}_{\star}$ relations is largely alleviated if similar methods are used to infer the galaxy properties. We find that relations relying on infrared wavelengths (e.g. 24 ${\rm \mu m}$, MIPS - 24, 70 and 160 ${\rm \mu m}$ or SPIRE - 250, 350, 500 ${\rm \mu m}$) have SFRs that exceed the intrinsic relation by 0.5 dex. Relations that rely on the spectral energy distribution fitting technique underpredict the SFRs at a fixed stellar mass by -0.5 dex at $z \sim 4$ but overpredict the measurements by 0.3 dex at $z \sim 1$. Relations relying on dust-corrected rest-frame UV luminosities, are flatter since they overpredict/underpredict SFRs for low/high star forming objects and yield deviations from the intrinsic relation from 0.10 dex to -0.13 dex at $z \sim 4$. We suggest that the severe tension between different observational studies can be broadly explained by the fact that different groups employ different techniques to infer their SFRs., Comment: 15 pages, 6 Figures, 4 tables, accepted for publication in MNRAS, Figures 1, 2, 3, 4 summarize results

The open source SKIRT Monte Carlo radiative transfer code has been used for more than 15 years to model the interaction between radiation and dust in various astrophysical systems. In this work, we present version 9 of the code, which has been substantially redesigned to support long-term objectives. We invite interested readers to participate in the development, testing and application of new features such as including gas media types in addition to dust, performing line transfer in addition to continuum radiation transfer, and modeling polarization by non-spherical dust grains aligned by magnetic fields. We describe the major challenges involved in preparing the code for these and other extensions, as well as their resolution, including a completely new treatment of wavelengths to support kinematics. SKIRT 9 properly runs over 400 handcrafted functional tests and successfully performs all relevant benchmarks. The source code and all documentation is publicly available for use and ready for further collaborative development., Comment: 24 pages, 14 figures, accepted for publication in Astronomy and Computing

We exploit EAGLE, a cosmological hydrodynamical simulation, to reproduce the selection of the observed sub-millimeter (submm) galaxy population by selecting the model galaxies at $z \geq 1$ with mock submm fluxes $S_{850} \geq 1$ mJy. There is a reasonable agreement between the galaxies within this sample and the properties of the observed submm population, such as their star formation rates (SFRs) at $z 1$ mJy. Thus, we investigate the nature of $z \geq 1$ highly star-forming Submm-Faint galaxies (i.e., $\dot M_* \geq 80$ Msol/yr but $S_{850}< 1$ mJy). We find they are similar to the model submm galaxies; being gas rich and hosting undermassive black holes, however they are typically lower mass ($M_* \sim 10^{10}$ Msol) and are at higher redshifts ($z>4$). These typically higher-$z$ galaxies show stronger evidence for having been triggered by major mergers, and critically, they are likely missed by current submm surveys due to their higher dust temperatures. This suggests a potentially even larger contribution to the SFR density at $z > 3$ from dust-obscured systems than implied by current observations., 15 pages, 9 figures, submitted to MNRAS

The EAGLE cosmological simulations reproduce the observed galaxy stellar mass function and many galaxy properties. In this work, we study the dust-related properties of present-day EAGLE galaxies through mock observations in the far-infrared and submm wavelength ranges obtained with the 3D dust radiative transfer code SKIRT. To prepare an EAGLE galaxy for radiative transfer processing, we derive a diffuse dust distribution from the gas particles and we re-sample the star-forming gas particles and the youngest star particles into star-forming regions that are assigned dedicated emission templates. We select a set of redshift-zero EAGLE galaxies that matches the K-band luminosity distribution of the galaxies in the Herschel Reference Survey (HRS), a volume-limited sample of about 300 normal galaxies in the Local Universe. We find overall agreement of the EAGLE dust scaling relations with those observed in the HRS, such as the dust-to-stellar mass ratio versus stellar mass and versus NUV-r colour relations. A discrepancy in the f_250/f_350 versus f_350/f_500 submm colour-colour relation implies that part of the simulated dust is insufficiently heated, likely because of limitations in our sub-grid model for star-forming regions. We also investigate the effect of adjusting the metal-to-dust ratio and the covering factor of the photodissociation regions surrounding the star-forming cores. We are able to constrain the important dust-related parameters in our method, informing the calculation of dust attenuation for EAGLE galaxies in the UV and optical domain., 19 pages, 16 figures, accepted for publication in MNRAS

We present the public data release of halo and galaxy catalogues extracted from the EAGLE suite of cosmological hydrodynamical simulations of galaxy formation. These simulations were performed with an enhanced version of the GADGET code that includes a modified hydrodynamics solver, time-step limiter and subgrid treatments of baryonic physics, such as stellar mass loss, element-by-element radiative cooling, star formation and feedback from star formation and black hole accretion. The simulation suite includes runs performed in volumes ranging from 25 to 100 comoving megaparsecs per side, with numerical resolution chosen to marginally resolve the Jeans mass of the gas at the star formation threshold. The free parameters of the subgrid models for feedback are calibrated to the redshift z=0 galaxy stellar mass function, galaxy sizes and black hole mass - stellar mass relation. The simulations have been shown to match a wide range of observations for present-day and higher-redshift galaxies. The raw particle data have been used to link galaxies across redshifts by creating merger trees. The indexing of the tree produces a simple way to connect a galaxy at one redshift to its progenitors at higher redshift and to identify its descendants at lower redshift. In this paper we present a relational database which we are making available for general use. A large number of properties of haloes and galaxies and their merger trees are stored in the database, including stellar masses, star formation rates, metallicities, photometric measurements and mock gri images. Complex queries can be created to explore the evolution of more than 10^5 galaxies, examples of which are provided in appendix. (abridged), Comment: 12 pages, 3 figures + appendices. The simulation data is available at http://www.eaglesim.org/database.php (comments welcome). Paper accepted to AstronomyComputing. Note there is a typo in the published version, Section 4.1, replace 3/2 with 2/3 for velocity dispersion

We present results of the detailed dust energy balance study for the seven large edge-on galaxies in the HEROES sample using 3D radiative transfer (RT) modelling. Based on available optical and near-infrared observations of the HEROES galaxies, we derive the 3D distribution of stars and dust in these galaxies. For the sake of uniformity, we apply the same technique to retrieve galaxy properties for the entire sample: we use a stellar model consisting of a S\'ersic bulge and three double-exponential discs (a superthin disc for a young stellar population and thin and thick discs for old populations). For the dust component, we adopt a double-exponential disc with the new THEMIS dust-grain model. We fit oligochromatic radiative transfer (RT) models to the optical and near-infrared images with the fitting algorithm FitSKIRT and do panchromatic simulations with the SKIRT code at wavelengths ranging from ultraviolet to submillimeter. We confirm the previously stated dust energy balance problem in galaxies: for the HEROES galaxies, the dust emission derived from our RT calculations underestimates the real observations by a factor 1.5-4 for all galaxies except NGC 973 and NGC 5907 (apparently, the latter galaxy has a more complex geometry than we used). The comparison between our RT simulations and the observations at mid-infrared-submillimeter wavelengths shows that most of our galaxies exhibit complex dust morphologies (possible spiral arms, star-forming regions, more extended dust structure in the radial and vertical directions). We suggest that, in agreement with the results from Saftly et al. (2015), the large- and small-scale structure is the most probable explanation for the dust energy balance problem., Comment: 53 pages, 31 figures, Accepted for publication in Astronomy & Astrophysics

The Monte Carlo method is the most popular technique to perform radiative transfer simulations in a general 3D geometry. The algorithms behind and acceleration techniques for Monte Carlo radiative transfer are discussed extensively in the literature, and many different Monte Carlo codes are publicly available. On the contrary, the design of a suite of components that can be used for the distribution of sources and sinks in radiative transfer codes has received very little attention. The availability of such models, with different degrees of complexity, has many benefits. For example, they can serve as toy models to test new physical ingredients, or as parameterised models for inverse radiative transfer fitting. For 3D Monte Carlo codes, this requires algorithms to efficiently generate random positions from 3D density distributions. We describe the design of a flexible suite of components for the Monte Carlo radiative transfer code SKIRT. The design is based on a combination of basic building blocks (which can be either analytical toy models or numerical models defined on grids or a set of particles) and the extensive use of decorators that combine and alter these building blocks to more complex structures. For a number of decorators, e.g. those that add spiral structure or clumpiness, we provide a detailed description of the algorithms that can be used to generate random positions. Advantages of this decorator-based design include code transparency, the avoidance of code duplication, and an increase in code maintainability. Moreover, since decorators can be chained without problems, very complex models can easily be constructed out of simple building blocks. Finally, based on a number of test simulations, we demonstrate that our design using customised random position generators is superior to a simpler design based on a generic black-box random position generator., Comment: 15 pages, 4 figures, accepted for publication in Astronomy and Computing

The radiation of stars heats dust grains in the diffuse interstellar medium and in star-forming regions in galaxies. Modelling this interaction provides information on dust in galaxies, a vital ingredient for their evolution. It is not straightforward to identify the stellar populations heating the dust, and to link attenuation to emission on a sub-galactic scale. Radiative transfer models are able to simulate this dust-starlight interaction in a realistic, three-dimensional setting. We investigate the dust heating mechanisms on a local and global galactic scale, using the Andromeda galaxy (M31) as our laboratory. We perform a series of panchromatic radiative transfer simulations of Andromeda with our code SKIRT. The high inclination angle of M31 complicates the 3D modelling and causes projection effects. However, the observed morphology and flux density are reproduced fairly well from UV to sub-millimeter wavelengths. Our model reveals a realistic attenuation curve, compatible with previous, observational estimates. We find that the dust in M31 is mainly (91 % of the absorbed luminosity) heated by the evolved stellar populations. The bright bulge produces a strong radiation field and induces non-local heating up to the main star-forming ring at 10 kpc. The relative contribution of unevolved stellar populations to the dust heating varies strongly with wavelength and with galactocentric distance.The dust heating fraction of unevolved stellar populations correlates strongly with NUV-r colour and specific star formation rate. These two related parameters are promising probes for the dust heating sources at a local scale., Comment: 16 pages + appendix, 13 figures, accepted for publication in A&A

The radiative transport of photons through arbitrary three-dimensional (3D) structures of dust is a challenging problem due to the anisotropic scattering of dust grains and strong coupling between different spatial regions. The radiative transfer problem in 3D is solved using Monte Carlo or Ray Tracing techniques as no full analytic solution exists for the true 3D structures. We provide the first 3D dust radiative transfer benchmark composed of a slab of dust with uniform density externally illuminated by a star. This simple 3D benchmark is explicitly formulated to provide tests of the different components of the radiative transfer problem including dust absorption, scattering, and emission. This benchmark includes models with a range of dust optical depths fully probing cases that are optically thin at all wavelengths to optically thick at most wavelengths. This benchmark includes solutions for the full dust emission including single photon (stochastic) heating as well as two simplifying approximations: One where all grains are considered in equilibrium with the radiation field and one where the emission is from a single effective grain with size-distribution-averaged properties. A total of six Monte Carlo codes and one Ray Tracing code provide solutions to this benchmark. Comparison of the results revealed that the global SEDs are consistent on average to a few percent for all but the scattered stellar flux at very high optical depths. The image results are consistent within 10%, again except for the stellar scattered flux at very high optical depths. The lack of agreement between different codes of the scattered flux at high optical depths is quantified for the first time. We provide the first 3D dust radiative transfer benchmark and validate the accuracy of this benchmark through comparisons between multiple independent codes and detailed convergence tests., Comment: 21 pages, 18 figures, Astronomy & Astrophysics, in press

We present mock optical images, broad-band and H$\alpha$ fluxes, and D4000 spectral indices for 30,145 galaxies from the EAGLE hydrodynamical simulation at redshift $z=0.1$, modelling dust with the SKIRT Monte Carlo radiative transfer code. The modelling includes a subgrid prescription for dusty star-forming regions, with both the subgrid obscuration of these regions and the fraction of metals in diffuse interstellar dust calibrated against far-infrared fluxes of local galaxies. The predicted optical colours as a function of stellar mass agree well with observation, with the SKIRT model showing marked improvement over a simple dust screen model. The orientation dependence of attenuation is weaker than observed because EAGLE galaxies are generally puffier than real galaxies, due to the pressure floor imposed on the interstellar medium. The mock H$\alpha$ luminosity function agrees reasonably well with the data, and we quantify the extent to which dust obscuration affects observed H$\alpha$ fluxes. The distribution of D4000 break values is bimodal, as observed. In the simulation, 20$\%$ of galaxies deemed `passive' for the SKIRT model, i.e. exhibiting D4000 $> 1.8$, are classified `active' when ISM dust attenuation is not included. The fraction of galaxies with stellar mass greater than $10^{10}$ M$_\odot$ that are deemed passive is slightly smaller than observed, which is due to low levels of residual star formation in these simulated galaxies. Colour images, fluxes and spectra of EAGLE galaxies are to be made available through the public EAGLE database., Comment: Accepted for publication in MNRAS, 33 pages, 21 figures, 3 Tables

We investigate the amount and spatial distribution of interstellar dust in edge-on spiral galaxies, using detailed radiative transfer modeling of a homogeneous sample of 12 galaxies selected from the CALIFA survey. Our automated fitting routine, FitSKIRT, was first validated against artificial data. This is done by simultaneously reproducing the SDSS $g$-, $r$-, $i$- and $z$-band observations of a toy model in order to combine the information present in the different bands. We show that this combined, oligochromatic fitting, has clear advantages over standard monochromatic fitting especially regarding constraints on the dust properties. We model all galaxies in our sample using a three-component model, consisting of a double exponential disc to describe the stellar and dust discs and using a S\'ersic profile to describe the central bulge. The full model contains 19 free parameters, and we are able to constrain all these parameters to a satisfactory level of accuracy without human intervention or strong boundary conditions. Apart from two galaxies, the entire sample can be accurately reproduced by our model. We find that the dust disc is about 75% more extended but only half as high as the stellar disc. The average face-on optical depth in the V-band is $0.76$ and the spread of $0.60$ within our sample is quite substantial, which indicates that some spiral galaxies are relatively opaque even when seen face-on., Comment: 18 pages, 6 figures, 4 tables, Accepted for publication in MNRAS

Context: The Monte Carlo method is the most widely used method to solve radiative transfer problems in astronomy, especially in a fully general 3D geometry. A crucial concept in any Monte Carlo radiative transfer code is the random generation of the next interaction location. In polarised Monte Carlo radiative transfer with aligned non-spherical grains, the nature of dichroism complicates the concept of optical depth. Aims: We investigate in detail the relation between optical depth and the optical properties and density of the attenuating medium in polarised Monte Carlo radiative transfer codes that take into account dichroic extinction. Methods: Based on solutions for the radiative transfer equation, we discuss the optical depth scale in polarised radiative transfer with spheroidal grains. We compare the dichroic optical depth to the extinction and total optical depth scale. Results: In a dichroic medium, the optical depth is not equal to the usual extinction optical depth, nor to the total optical depth. For representative values of the optical properties of dust grains, the dichroic optical depth can differ from the extinction or total optical depth by several ten percent. A closed expression for the dichroic optical depth cannot be given, but it can be derived efficiently through an algorithm that is based on the analytical result corresponding to elongated grains with a uniform grain alignment. Conclusions: Optical depth is more complex in dichroic media than in systems without dichroic attenuation, and this complexity needs to be considered when generating random free path lengths in Monte Carlo radiative transfer simulations. There is no benefit in using approximations instead of the dichroic optical depth., 9 pages, 3 figures, accepted for publication in A&A

Biasing or importance sampling is a powerful technique in Monte Carlo radiative transfer, and can be applied in different forms to increase the accuracy and efficiency of simulations. One of the drawbacks of the use of biasing is the potential introduction of large weight factors. We discuss a general strategy, composite biasing, to suppress the appearance of large weight factors. We use this composite biasing approach for two different problems faced by current state-of-the-art Monte Carlo radiative transfer codes: the generation of photon packages from multiple components, and the penetration of radiation through high optical depth barriers. In both cases, the implementation of the relevant algorithms is trivial and does not interfere with any other optimisation techniques. Through simple test models, we demonstrate the general applicability, accuracy and efficiency of the composite biasing approach. In particular, for the penetration of high optical depths, the gain in efficiency is spectacular for the specific problems that we consider: in simulations with composite path length stretching, high accuracy results are obtained even for simulations with modest numbers of photon packages, while simulations without biasing cannot reach convergence, even with a huge number of photon packages., 12 pages, accepted for publication in A&A

Interstellar dust in galaxies can be traced either through its extinction effects on the star light, or through its thermal emission at infrared wavelengths. Recent radiative transfer studies of several nearby edge-on galaxies have found an apparent inconsistency in the dust energy balance: the radiative transfer models that successfully explain the optical extinction underestimate the observed fluxes by an average factor of three. We investigate the dust energy balance for IC4225 and NGC5166, two edge-on spiral galaxies observed by the Herschel Space Observatory in the frame of the H-ATLAS survey. We start from models which were constrained from optical data and extend them to construct the entire spectral energy distribution of our galaxies. These predicted values are subsequently compared to the observed far-infrared fluxes. We find that including a young stellar population in the modelling is necessary as it plays a non-negligible part in the heating of the dust grains. While the modelling approach for both galaxies is nearly identical, we find two very different results. As is often seen in other edge-on spiral galaxies, the far-infrared emission of our radiative transfer model of IC4225 underestimates the observed fluxes by a factor of about three. For NGC5166 on the other hand, we find that both the predicted spectral energy distribution as well as the simulated images match the observations particularly well. We explore possible reasons for this difference and conclude that it is unlikely that one single mechanism is the cause of the dust energy balance problem in spiral galaxies. We discuss the different approaches that can be considered in order to get a conclusive answer on the origin this discrepancy., 12 pages, 4 figures, accepted for publication in MNRAS

In light of the recent detection of direct evidence for the formation of Kelvin-Helmholtz instabilities in the Orion nebula, we expand upon previous modelling efforts by numerically simulating the shear-flow driven gas and dust dynamics in locations where the H$_{II}$ region and the molecular cloud interact. We aim to directly confront the simulation results with the infrared observations. Methods: To numerically model the onset and full nonlinear development of the Kelvin-Helmholtz instability we take the setup proposed to interpret the observations, and adjust it to a full 3D hydrodynamical simulation that includes the dynamics of gas as well as dust. A dust grain distribution with sizes between 5-250 nm is used, exploiting the gas+dust module of the MPI-AMRVAC code, in which the dust species are represented by several pressureless dust fluids. The evolution of the model is followed well into the nonlinear phase. The output of these simulations is then used as input for the SKIRT dust radiative transfer code to obtain infrared images at several stages of the evolution, which can be compared to the observations. Results: We confirm that a 3D Kelvin-Helmholtz instability is able to develop in the proposed setup, and that the formation of the instability is not inhibited by the addition of dust. Kelvin-Helmholtz billows form at the end of the linear phase, and synthetic observations of the billows show striking similarities to the infrared observations. It is pointed out that the high density dust regions preferentially collect on the flanks of the billows. To get agreement with the observed Kelvin-Helmholtz ripples, the assumed geometry between the background radiation, the billows and the observer is seen to be of critical importance., 8 pages, 10 figures

The interstellar dust content in galaxies can be traced in extinction at optical wavelengths, or in emission in the far-infrared. Several studies have found that radiative transfer models that successfully explain the optical extinction in edge-on spiral galaxies generally underestimate the observed FIR/submm fluxes by a factor of about three. In order to investigate this so-called dust energy balance problem, we use two Milky Way-like galaxies produced by high-resolution hydrodynamical simulations. We create mock optical edge-on views of these simulated galaxies (using the radiative transfer code SKIRT), and we then fit the parameters of a basic spiral galaxy model to these images (using the fitting code FitSKIRT). The basic model includes smooth axisymmetric distributions along a S\'ersic bulge and exponential disc for the stars, and a second exponential disc for the dust. We find that the dust mass recovered by the fitted models is about three times smaller than the known dust mass of the hydrodynamical input models. This factor is in agreement with previous energy balance studies of real edge-on spiral galaxies. On the other hand, fitting the same basic model to less complex input models (e.g. a smooth exponential disc with a spiral perturbation or with random clumps), does recover the dust mass of the input model almost perfectly. Thus it seems that the complex asymmetries and the inhomogeneous structure of real and hydrodynamically simulated galaxies are a lot more efficient at hiding dust than the rather contrived geometries in typical quasi-analytical models. This effect may help explain the discrepancy between the dust emission predicted by radiative transfer models and the observed emission in energy balance studies for edge-on spiral galaxies., Comment: 9 pages, 5 figures, accepted for publication in A&A

Polarization is an important tool to further the understanding of interstellar dust and the sources behind it. In this paper we describe our implementation of polarization that is due to scattering of light by spherical grains and electrons in the dust Monte Carlo radiative transfer code SKIRT. In contrast to the implementations of other Monte Carlo radiative transfer codes, ours uses co-moving reference frames that rely solely on the scattering processes. It fully supports the peel-off mechanism that is crucial for the efficient calculation of images in 3D Monte Carlo codes. We develop reproducible test cases that push the limits of our code. The results of our program are validated by comparison with analytically calculated solutions. Additionally, we compare results of our code to previously published results. We apply our method to models of dusty spiral galaxies at near-infrared and optical wavelengths. We calculate polarization degree maps and show them to contain signatures that trace characteristics of the dust arms independent of the inclination or rotation of the galaxy., 15 pages, 10 figures, accepted by Astronomy and Astrophysics

Context: Dust reprocesses about half of the stellar radiation in galaxies. The thermal re-emission by dust of absorbed energy is considered driven merely by young stars and, consequently, often applied to trace the star formation rate in galaxies. Recent studies have argued that the old stellar population might anticipate a non-negligible fraction of the radiative dust heating. Aims: In this work, we aim to analyze the contribution of young (< 100 Myr) and old (~ 10 Gyr) stellar populations to radiative dust heating processes in the nearby grand-design spiral galaxy M51 using radiative transfer modeling. High-resolution 3D radiative transfer (RT) models are required to describe the complex morphologies of asymmetric spiral arms and clumpy star-forming regions and model the propagation of light through a dusty medium. Methods: In this paper, we present a new technique developed to model the radiative transfer effects in nearby face-on galaxies. We construct a high-resolution 3D radiative transfer model with the Monte-Carlo code SKIRT accounting for the absorption, scattering and non-local thermal equilibrium (NLTE) emission of dust in M51. The 3D distribution of stars is derived from the 2D morphology observed in the IRAC 3.6 {\mu}m, GALEX FUV, H{\alpha} and MIPS 24 {\mu}m wavebands, assuming an exponential vertical distribution with an appropriate scale height. The dust geometry is constrained through the far-ultraviolet (FUV) attenuation, which is derived from the observed total-infrared-to-far-ultraviolet luminosity ratio. The stellar luminosity, star formation rate and dust mass have been scaled to reproduce the observed stellar spectral energy distribution (SED), FUV attenuation and infrared SED. (abridged), Comment: Accepted for publication in Astronomy & Astrophysics, 23 pages, 15 figures

We investigate the dust energy balance for the edge-on galaxy IC 2531, one of the seven galaxies in the HEROES sample. We perform a state-of-the-art radiative transfer modelling based, for the first time, on a set of optical and near-infrared galaxy images. We show that taking into account near-infrared imaging in the modelling significantly improves the constraints on the retrieved parameters of the dust content. We confirm the result from previous studies that including a young stellar population in the modelling is important for explaining the observed stellar energy distribution. However, the discrepancy between the observed and modelled thermal emission at far-infrared wavelengths, the so-called dust energy balance problem, is still present: the model underestimates the observed fluxes by a factor of about two. We compare two different dust models, and find that dust parameters and thus the spectral energy distribution in the infrared domain are sensitive to the adopted dust model. In general, the THEMIS model reproduces the observed emission in the infrared wavelength domain better than the popular Zubko et al. BARE-GR-S model. Our study of IC 2531 is a pilot case for detailed and uniform radiative transfer modelling of the entire HEROES sample, which will shed more light on the strength and origins of the dust energy balance problem., Comment: 15 pages, 8 figures, Accepted for publication in Astronomy & Astrophysics

We introduce the Virgo Consortium's EAGLE project, a suite of hydrodynamical simulations that follow the formation of galaxies and black holes in representative volumes. We discuss the limitations of such simulations in light of their finite resolution and poorly constrained subgrid physics, and how these affect their predictive power. One major improvement is our treatment of feedback from massive stars and AGN in which thermal energy is injected into the gas without the need to turn off cooling or hydrodynamical forces, allowing winds to develop without predetermined speed or mass loading factors. Because the feedback efficiencies cannot be predicted from first principles, we calibrate them to the z~0 galaxy stellar mass function and the amplitude of the galaxy-central black hole mass relation, also taking galaxy sizes into account. The observed galaxy mass function is reproduced to $\lesssim 0.2$ dex over the full mass range, $10^8 < M_*/M_\odot \lesssim 10^{11}$, a level of agreement close to that attained by semi-analytic models, and unprecedented for hydrodynamical simulations. We compare our results to a representative set of low-redshift observables not considered in the calibration, and find good agreement with the observed galaxy specific star formation rates, passive fractions, Tully-Fisher relation, total stellar luminosities of galaxy clusters, and column density distributions of intergalactic CIV and OVI. While the mass-metallicity relations for gas and stars are consistent with observations for $M_* \gtrsim 10^9 M_\odot$, they are insufficiently steep at lower masses. The gas fractions and temperatures are too high for clusters of galaxies, but for groups these discrepancies can be resolved by adopting a higher heating temperature in the subgrid prescription for AGN feedback. EAGLE constitutes a valuable new resource for studies of galaxy formation., Accepted for publication in MNRAS. V2: Minor changes. For images and videos, see http://eagle.strw.leidenuniv.nl/ and http://icc.dur.ac.uk/Eagle/

We present redshift-zero synthetic dust-aware observations for the 45 Milky Way-mass simulated galaxies of the ARTEMIS project, calculated with the SKIRT radiative transfer code. The post-processing procedure includes components for star-forming regions, stellar sources, and diffuse dust. We produce and publicly release realistic high-resolution images for 50 commonly-used broadband filters from ultraviolet to sub-millimetre wavelengths and for 18 different viewing angles. We compare the simulated ARTEMIS galaxies to observed galaxies in the DustPedia database with similar stellar mass and star formation rate, and to synthetic observations of the simulated galaxies of the Auriga project produced in previous work using a similar post-processing technique. In all cases, global galaxy properties are derived using SED fitting. We find that, similar to Auriga, the post-processed ARTEMIS galaxies generally reproduce the observed scaling relations for global fluxes and physical properties, although dust extinction at FUV/UV wavelengths is underestimated and representative dust temperatures are lower than observed. At a resolved scale, we compare multi-wavelength non-parametric morphological properties of selected disc galaxies across the data sets. We find that the ARTEMIS galaxies largely reproduce the observed morphological trends as a function of wavelength, although they appear to be more clumpy and less symmetrical than observed. We note that the ARTEMIS and Auriga galaxies occupy adjacent regions in the specific star formation versus stellar mass plane, so that the synthetic observation data sets supplement each other., Accepted for publication by MNRAS

We present dust-attenuated and dust emission fluxes for sufficiently resolved galaxies in the EAGLE suite of cosmological hydrodynamical simulations, calculated with the SKIRT radiative transfer code. The post-processing procedure includes specific components for star formation regions, stellar sources, and diffuse dust, and takes into account stochastic heating of dust grains to obtain realistic broad-band fluxes in the wavelength range from ultraviolet to sub-millimeter. The mock survey includes nearly half a million simulated galaxies with stellar masses above 10^8.5 solar masses across six EAGLE models. About two thirds of these galaxies, residing in 23 redshift bins up to z=6, have a sufficiently resolved metallic gas distribution to derive meaningful dust attenuation and emission, with the important caveat that the same dust properties were used at all redshifts. These newly released data complement the already publicly available information about the EAGLE galaxies, which includes intrinsic properties derived by aggregating the properties of the smoothed particles representing matter in the simulation. We further provide an open source framework of Python procedures for post-processing simulated galaxies with the radiative transfer code SKIRT. The framework allows any third party to calculate synthetic images, SEDs, and broadband fluxes for EAGLE galaxies, taking into account the effects of dust attenuation and emission., 20 pages, 10 figures, accepted for publication in ApJS