Submillimetre galaxies (SMGs) have long posed a challenge for theorists, and self-consistently reproducing the properties of the SMG population in a large-volume cosmological hydrodynamical simulation has not yet been achieved. In this work, we employ a relatively simple method based on previous work to predict the submm flux densities of simulated SMGs drawn from cosmological simulations from the Illustris and IllustrisTNG projects and compare the predicted number counts with observations. We find that the predicted SMG number counts based on IllustrisTNG are significantly less than observed (more than 1 dex at $S_{850} \gtrsim 4$ mJy). The simulation from the original Illustris project yields more SMGs than IllustrisTNG: the predicted counts are consistent with those observed at both $S_{850} \lesssim 5$ mJy and $S_{850} \gtrsim 9$ mJy and only a factor of $\sim 2$ lower than observed at intermediate flux densities. We demonstrate that IllustrisTNG hosts fewer SMGs than Illustris because in the former, high-mass ($M_{\star} \sim 10^{11} \, {\rm M}_{\odot}$) $z \sim 2-3$ galaxies have lower dust masses and star formation rates (SFRs) than in Illustris owing to differences in the sub-grid models for stellar or/and active galactic nucleus (AGN) feedback between the two simulations (we unfortunately cannot isolate the specific cause(s) post hoc). Our results demonstrate that because our method enables predicting SMG number counts in post-processing with a negligible computational expense, SMGs can provide useful constraints for tuning sub-grid models in future large-volume cosmological simulations. Higher resolution, which would lead to stronger SFR enhancements in starbursts, could at least partially reconcile the discrepancy between the IllustrisTNG SMG number counts and those observed.

We compare predictions for galaxy-galaxy lensing profiles and clustering from the Henriques et al. (2015) public version of the Munich semi-analytical model of galaxy formation (SAM) and the IllustrisTNG suite, primarily TNG300, with observations from KiDS+GAMA and SDSS-DR7 using four different selection functions for the lenses (stellar mass, stellar mass and group membership, stellar mass and isolation criteria, stellar mass and colour). We find that this version of the SAM does not agree well with the current data for stellar mass-only lenses with $M_\ast > 10^{11}\,M_\odot$. By decreasing the merger time for satellite galaxies as well as reducing the radio-mode AGN accretion efficiency in the SAM, we obtain better agreement, both for the lensing and the clustering, at the high mass end. We show that the new model is consistent with the signals for central galaxies presented in Velliscig et al. (2017). Turning to the hydrodynamical simulation, TNG300 produces good lensing predictions, both for stellar mass-only ($\chi^2 = 1.81$ compared to $\chi^2 = 7.79$ for the SAM), and locally brightest galaxies samples ($\chi^2 = 3.80$ compared to $\chi^2 = 5.01$). With added dust corrections to the colours it matches the SDSS clustering signal well for red low mass galaxies. We find that both the SAMs and TNG300 predict $\sim 50\,\%$ excessive lensing signals for intermediate mass red galaxies with $10.2 < \log_{10} M_\ast [ M_\odot ] < 11.2$ at $r \approx 0.6\,h^{-1}\,\mathrm{Mpc}$, which require further theoretical development.

Hydrodynamical simulations of galaxy formation and evolution attempt to fully model the physics that shapes galaxies. The agreement between the morphology of simulated and real galaxies, and the way the morphological types are distributed across galaxy scaling relations are important probes of our knowledge of galaxy formation physics. Here we propose an unsupervised deep learning approach to perform a stringent test of the fine morphological structure of galaxies coming from the Illustris and IllustrisTNG (TNG100 and TNG50) simulations against observations from a subsample of the Sloan Digital Sky Survey. Our framework is based on PixelCNN, an autoregressive model for image generation with an explicit likelihood. We adopt a strategy that combines the output of two PixelCNN networks in a metric that isolates the fine morphological details of galaxies from the sky background. We are able to \emph{quantitatively} identify the improvements of IllustrisTNG, particularly in the high-resolution TNG50 run, over the original Illustris. However, we find that the fine details of galaxy structure are still different between observed and simulated galaxies. This difference is driven by small, more spheroidal, and quenched galaxies which are globally less accurate regardless of resolution and which have experienced little improvement between the three simulations explored. We speculate that this disagreement, that is less severe for quenched disky galaxies, may stem from a still too coarse numerical resolution, which struggles to properly capture the inner, dense regions of quenched spheroidal galaxies.

Dynamically relaxed galaxy clusters have long played a role in galaxy cluster studies because it is thought their properties can be reconstructed more precisely and with less systematics. As relaxed clusters are desirable, there exist a plethora of criteria for classifying a galaxy cluster as relaxed. In this work, we examine $9$ commonly used observational and theoretical morphological metrics extracted from $54,000$ Mock-X synthetic X-ray images of galaxy clusters taken from the IllustrisTNG, BAHAMAS and MACSIS simulation suites. We find that the simulated criteria distributions are in reasonable agreement with the observed distributions. Many criteria distributions evolve as a function of redshift, cluster mass, numerical resolution and subgrid physics, limiting the effectiveness of a single relaxation threshold value. All criteria are positively correlated with each other, however, the strength of the correlation is sensitive to redshift, mass and numerical choices. Driven by the intrinsic scatter inherent to all morphological metrics and the arbitrary nature of relaxation threshold values, we find the consistency of relaxed subsets defined by the different metrics to be relatively poor. Therefore, the use of relaxed cluster subsets introduces significant selection effects that are non-trivial to resolve.

The past decade has seen significant progress in understanding galaxy formation and evolution using large-scale cosmological simulations. While these simulations produce galaxies in overall good agreement with observations, they employ different sub-grid models for galaxies and supermassive black holes (BHs). We investigate the impact of the sub-grid models on the BH mass properties of the Illustris, TNG100, TNG300, Horizon-AGN, EAGLE, and SIMBA simulations, focusing on the M_BH-M_star relation and the BH mass function. All simulations predict tight M_BH-M_star relations, and struggle to produce the lowest (M_BH< 10^7.5 Msun) in galaxies of M_star~10^10.5-10^11.5 Msun. While the time evolution of the mean M_BH-M_star relation is mild (<1 dex in BH mass for 0<z<5) for all the simulations, its linearity (shape) and normalization varies from simulation to simulation. The strength of SN feedback has a large impact on the linearity and time evolution for M_star<10^10.5 Msun. We find that the low-mass end is a good discriminant of the simulation models, and highlights the need for new observational constraints. At the high-mass end, strong AGN feedback can suppress the time evolution of the relation normalization. Compared with the observations of the local universe, we find an excess of BHs with M_BH>10^9 Msun in most of the simulations. The BH mass function is dominated by efficiently accreting BHs (log10 f_Edd >-2$) at high redshifts, and transitions progressively from the high-mass to the low-mass end to be governed by inactive BHs. The transition time and the contribution of active BHs are different among the simulations, and can be used to evaluate models against observations.

We measure power spectrum response functions in the presence of baryonic physical processes using separate universe simulations with the IllustrisTNG galaxy formation model. The response functions describe how the small-scale power spectrum reacts to long-wavelength perturbations and they can be efficiently measured with the separate universe technique by absorbing the effects of the long modes into a modified cosmology. Specifically, we focus on the total first-order matter power spectrum response to an isotropic density fluctuation $R_1(k,z)$, which is fully determined by the logarithmic derivative of the nonlinear matter power spectrum ${\rm dln}P_m(k,z)/{\rm dln}k$ and the growth-only response function $G_1(k,z)$. We find that $G_1(k,z)$ is not affected by the baryonic physical processes in the simulations at redshifts $z < 3$ and on all scales probed ($k \lesssim 15h/{\rm Mpc}$, i.e. length scales $\gtrsim 0.4 {\rm Mpc}/h$). In practice, this implies that the power spectrum fully specifies the baryonic dependence of its response function. Assuming an idealized lensing survey setup, we evaluate numerically the baryonic impact on the squeezed-lensing bispectrum and the lensing super-sample power spectrum covariance, which are given in terms of responses. Our results show that these higher-order lensing statistics can display varying levels of sensitivity to baryonic effects compared to the power spectrum, with the squeezed-bispectrum being the least sensitive. We also show that ignoring baryonic effects on lensing covariances slightly overestimates the error budget (and is therefore conservative from the point of view of parameter error bars) and likely has negligible impact on parameter biases in inference analyses.

We perform a series of controlled N-body simulations to study realizations of the recently discovered Antlia 2 galaxy in cold dark matter (CDM) and self-interacting dark matter (SIDM) scenarios. Our simulations contain six benchmark models, where we vary the initial halo concentration and the self-scattering cross section. We adopt well-motivated initial stellar and halo masses, and our fiducial orbit has a small pericenter. After evolving in the Milky Way's tidal field, the simulated galaxies experience significant mass loss and their stellar distributions expand accordingly. These tidal effects are more prominent if the initial halo concentration is lower and if the self-scattering cross section is larger. Our results show that Antlia 2-like galaxies could be realized in CDM if the halo concentration is low and the stellar distribution is diffuse at the infall time, while these conditions could be relaxed in SIDM. We also find all the simulated galaxies predict approximately the same stellar velocity dispersion after imposing selection criteria for stellar particles. This has important implications for testing dark matter models using tidally disturbed systems.

We propose two effective parameters that fully characterise galactic-scale structure formation at high redshifts ($z\gtrsim5$) for a variety of dark matter (DM) models that have a primordial cutoff in the matter power spectrum. Our description is within the recently proposed ETHOS framework and includes standard thermal Warm DM (WDM) and models with dark acoustic oscillations (DAOs). To define and explore this parameter space, we use high-redshift zoom-in simulations that cover a wide range of non-linear scales from those where DM should behave as CDM ($k\sim10\,h\,{\rm Mpc}^{-1}$), down to those characterised by the onset of galaxy formation ($k\sim500\,h\,{\rm Mpc}^{-1}$). We show that the two physically motivated parameters $h_{\rm peak}$ and $k_{\rm peak}$, the amplitude and scale of the first DAO peak, respectively, are sufficient to parametrize the linear matter power spectrum and classify the DM models as belonging to effective non-linear structure formation regions. These are defined by their relative departure from Cold DM ($k_{\rm peak}\rightarrow\infty$) and WDM ($h_{\rm peak}=0$) according to the non-linear matter power spectrum and halo mass function. We identify a region where the DAOs still leave a distinct signature from WDM down to $z=5$, while a large part of the DAO parameter space is shown to be degenerate with WDM. Our framework can then be used to seamlessly connect a broad class of particle DM models to their structure formation properties at high redshift without the need of additional $N$-body simulations.

We perform a suite of hydrodynamic simulations to investigate how initial density profiles of giant molecular clouds (GMCs) affect their subsequent evolution. We find that the star formation duration and integrated star formation efficiency of the whole clouds are not sensitive to the choice of different profiles but are mainly controlled by the interplay between gravitational collapse and stellar feedback. Despite this similarity, GMCs with different profiles show dramatically different modes of star formation. For shallower profiles, GMCs first fragment into many self-gravitation cores and form sub-clusters that distributed throughout the entire clouds. These sub-clusters are later assembled ``hierarchically'' to central clusters. In contrast, for steeper profiles, a massive cluster is quickly formed at the center of the cloud and then gradually grows its mass via gas accretion. Consequently, central clusters that emerged from clouds with shallower profiles are less massive and show less rotation than those with the steeper profiles. This is because 1) a significant fraction of mass and angular momentum in shallower profiles is stored in the orbital motion of the sub-clusters that are not able to merge into the central clusters 2) frequent hierarchical mergers in the shallower profiles lead to further losses of mass and angular momentum via violent relaxation and tidal disruption. Encouragingly, the degree of cluster rotations in steeper profiles is consistent with recent observations of young and intermediate-age clusters. We speculate that rotating globular clusters are likely formed via an ``accretion'' mode from centrally-concentrated clouds in the early Universe.

We present predictions for high redshift ($z=2-10$) galaxy populations based on the IllustrisTNG simulation suite and a full Monte Carlo dust radiative transfer post-processing. Specifically, we discuss the ${\rm H}_{\alpha}$ and ${\rm H}_{\beta}$ + $[\rm O \,III]$ luminosity functions up to $z=8$. The predicted ${\rm H}_{\beta}$ + $[\rm O \,III]$ luminosity functions are consistent with present observations at $z\lesssim 3$ with $\lesssim 0.1\,{\rm dex}$ differences in luminosities. However, the predicted ${\rm H}_{\alpha}$ luminosity function is $\sim 0.3\,{\rm dex}$ dimmer than the observed one at $z\simeq 2$. Furthermore, we explore continuum spectral indices, the Balmer break at $4000$\AA (D4000) and the UV continuum slope $\beta$. The median D4000 versus sSFR relation predicted at $z=2$ is in agreement with the local calibration despite a different distribution pattern of galaxies in this plane. In addition, we reproduce the observed $A_{\rm UV}$ versus $\beta$ relation and explore its dependence on galaxy stellar mass, providing an explanation for the observed complexity of this relation. We also find a deficiency in heavily attenuated, UV red galaxies in the simulations. Finally, we provide predictions for the dust attenuation curves of galaxies at $z=2-6$ and investigate their dependence on galaxy colors and stellar masses. The attenuation curves are steeper in galaxies at higher redshifts, with bluer colors, or with lower stellar masses. We attribute these predicted trends to dust geometry. Overall, our results are consistent with present observations of high redshift galaxies. Future JWST observations will further test these predictions.

We use the high-resolution TNG50 cosmological magnetohydrodynamical simulation to explore the properties and origin of cold circumgalactic medium (CGM) gas around massive galaxies (M* > 10^11 Msun) at intermediate redshift (z~0.5). We discover a significant abundance of small-scale, cold gas structure in the CGM of 'red and dead' elliptical systems, as traced by neutral HI and MgII. Halos can host tens of thousands of discrete absorbing cloudlets, with sizes of order a kpc or smaller. With a Lagrangian tracer analysis, we show that cold clouds form due to strong drho/rho >> 1 gas density perturbations which stimulate thermal instability. These local overdensities trigger rapid cooling from the hot virialized background medium at ~10^7 K to radiatively inefficient ~10^4 K clouds, which act as cosmologically long-lived, 'stimulated cooling' seeds in a regime where the global halo does not satisfy the classic tcool/tff < 10 criterion. Furthermore, these small clouds are dominated by magnetic rather than thermal pressure, with plasma beta << 1, suggesting that magnetic fields may play an important role. The number and total mass of cold clouds both increase with resolution, and the ~8x10^4 Msun cell mass of TNG50 enables the ~few hundred pc, small-scale CGM structure we observe to form. Finally, we contrast against observations from the COS-LRG, LRG-RDR, COS-Halos, and SDSS LRG surveys. We broadly find that our recent, high-resolution cosmological simulations produce sufficiently high covering fractions of extended, cold gas as observed to surround massive galaxies.

Low surface brightness galactic stellar halos provide a challenging but promising path towards unraveling the past assembly histories of individual galaxies. Here, we present detailed comparisons between the stellar halos of Milky Way-mass disk galaxies observed as part of the Dragonfly Nearby Galaxies Survey (DNGS) and stellar mass-matched galaxies in the TNG100 run of the IllustrisTNG project. We produce stellar mass maps as well as mock $g$ and $r$-band images for randomly oriented simulated galaxies, convolving the latter with the Dragonfly PSF and taking care to match the background noise, surface brightness limits and spatial resolution of DNGS. We measure azimuthally averaged stellar mass density and surface brightness profiles, and find that the DNGS galaxies generally have less stellar mass (or light) at large radii (>20 kpc) compared to their mass-matched TNG100 counterparts, and that simulated galaxies with similar surface density profiles tend to have low accreted mass fractions for their stellar mass. We explore potential solutions to this apparent "missing outskirts problem" by implementing several ad-hoc adjustments within TNG100 at the stellar particle level. Although we are unable to identify any single adjustment that fully reconciles the differences between the observed and simulated galaxy outskirts, we find that artificially delaying the disruption of satellite galaxies and reducing the spatial extent of in-situ stellar populations result in improved matches between the outer profile shapes and stellar halo masses, respectively. Further insight can be achieved with higher resolution simulations that are able to better resolve satellite accretion, and with larger samples of observed galaxies.

Supermassive black holes (SMBHs) which reside at the centres of galaxies can inject vast amounts of energy into the surrounding gas and are thought to be a viable mechanism to quench star-formation in massive galaxies. Here we study the $10^{9-12.5}\,\mathrm{M_\odot}$ stellar mass central galaxy population of the IllustrisTNG simulation, specifically the TNG100 and TNG300 volumes at $z=0$, and show how the three components - SMBH, galaxy, and circumgalactic medium (CGM) - are interconnected in their evolution. We find that gas entropy is a sensitive diagnostic of feedback injection. In particular, we demonstrate how the onset of the low-accretion BH feedback mode, realised in the IllustrisTNG model as a kinetic, BH-driven wind, leads not only to star-formation quenching at stellar masses $\gtrsim10^{10.5}\,\mathrm{M_\odot}$ but also to a change in thermodynamic properties of the (non-star-forming) gas, both within the galaxy and beyond. The IllustrisTNG kinetic feedback from SMBHs increases the average gas entropy, within the galaxy and in the CGM, lengthening typical gas cooling times from $10-100\, \mathrm{Myr}$ to $1-10\,\mathrm{Gyr}$, effectively ceasing ongoing star-formation and inhibiting radiative cooling and future gas accretion. In practice, the same AGN feedback channel is simultaneously `ejective' and `preventative' and leaves an imprint on the temperature, density, entropy, and cooling times also in the outer reaches of the gas halo, up to distances of several hundred kiloparsecs. In the IllustrisTNG model, a long-lasting quenching state can occur for a heterogeneous CGM, whereby the hot and dilute CGM gas of quiescent galaxies contains regions of low-entropy gas with short cooling times.

We investigate the influence of environment on the cold-gas properties of galaxies at z=0 within the TNG100 cosmological, magnetohydrodynamic simulation, part of the IllustrisTNG suite. We extend previous post-processing methods for breaking gas cells into their atomic and molecular phases, and build detailed mocks to comprehensively compare to the latest surveys of atomic hydrogen (HI) in nearby galaxies, namely ALFALFA and xGASS. We use TNG100 to explore the HI content, star formation activity, and angular momentum of satellite galaxies, each as a function of environment, and find that satellites are typically a factor of ~3 poorer in HI than centrals of the same stellar mass, with the exact offset depending sensitively on parent halo mass. Due to the large physical scales on which HI measurements are made (~45--245 kpc), contributions from gas not bound to the galaxy of interest but in the same line of sight crucially lead to larger HI mass measurements in the mocks in many cases, ultimately aligning with observations. This effect is mass-dependent and naturally greater for satellites than centrals, as satellites are never isolated by definition. We also show that HI stripping in TNG100 satellites is closely accompanied by quenching, in tension with observational data that instead favour that HI is preferentially stripped before star formation is reduced.

We study the stellar morphological evolution of disk galaxies within clusters in the TNG50 and TNG100 runs from the IllustrisTNG simulation suite. We select satellites of masses $10^{9.7}<=M_{*,z=0}/M_{sun}<=10^{11.6}$ residing in clusters of total masses $10^{14}<=M_{\text{200c,z=0}}/M_{sun}<10^{14.6}$ at z=0 and study those that were disks at accretion according to a kinematic morphology indicator (the circularity fraction). The galaxies' histories are traced from the time of accretion to $z=0$ and compared to a control sample of central galaxies mass-matched at the time of accretion. Most cluster disks become non-disky by z=0, in stark contrast with the control disks, of which a significant fraction remains disky over the same timescales. The transformation to non-disky morphologies is accompanied by gas removal and star formation quenching for both cluster and control galaxies. However, cluster disks that become non-disky by z=0 have lost dark matter (DM) mass and show little growth or a loss of stellar mass, whereas the corresponding control disks show significant growth in both components. Most cluster satellites change their morphologies on similar timescales regardless of stellar mass, in ~0.5-4 Gyr after accretion. Cluster disks that have had more numerous and closer pericentric passages show the largest change in morphology. Morphological change in both cluster and control disks requires the presence of a gravitational perturbation to drive stellar orbits to non-disky configurations, along with gas removal/heating to prevent replenishment of the disk through continued star-formation. For cluster disks, the perturbation is in the form of impulsive tidal shocking at pericentres and not tidal stripping of the outer disk stellar material, whereas for control disks, a combination of mergers and AGN feedback appears to be the key driving force behind morphological transformations.

We study evolution of self-interacting dark matter (SIDM) subhalos in the Milky Way (MW) tidal field. The interaction between the subhalos and the MW's tides lead to more diverse dark matter distribution in the inner region, compared to their cold dark matter counterparts. We test this scenario with two MW satellite galaxies, Draco and Fornax, opposite extremes in the inner dark matter content, and find that they can be accommodated within the SIDM model proposed to explain the diverse rotation curves of spiral galaxies in the field.

Satellites constitute an important fraction of the overall galaxy population and are believed to form in dark matter subhalos. Here we use the cosmological hydrodynamic simulation TNG100 to investigate how the formation histories of subhalos affect the properties and evolution of their host galaxies. We use a scaled formation time ($a_{\rm nf}$) to characterize the mass assembly histories of the subhalos before they are accreted by massive host halos. We find that satellite galaxies in young subhalos (low $a_{\rm nf}$) are less massive and more gas rich, and have stronger star formation and a higher fraction of ex situ stellar mass than satellites in old subhalos (high $a_{\rm nf}$). Furthermore, these low $a_{\rm nf}$ satellites require longer timescales to be quenched as a population than the high $a_{\rm nf}$ counterparts. We find very different merger histories between satellites in fast accretion (FA, $a_{\rm nf}<1.3$) and slow accretion (SA, $a_{\rm nf}>1.3$) subhalos. For FA satellites, the galaxy merger frequency dramatically increases just after accretion, which enhances the star formation at accretion. While, for SA satellites, the mergers occur smoothly and continuously across the accretion time. Moreover, mergers with FA satellites happen mainly after accretion, while a contrary trend is found for SA satellites. Our results provide insight into the evolution and star formation quenching of the satellite population.

Hot gaseous atmospheres that permeate galaxies and extend far beyond their stellar distribution, where they are commonly referred to as the circumgalactic medium (CGM), imprint important information about feedback processes powered by the stellar populations of galaxies and their central supermassive black holes (SMBH). In this work we study the properties of this hot X-ray emitting medium using the IllustrisTNG cosmological simulations. We analyse their mock X-ray spectra, obtained from the diffuse and metal-enriched gas in TNG100 and TNG50, and compare the results with X-ray observations of nearby early-type galaxies. The simulations reproduce the observed X-ray luminosities ($L_{\rm X}$) and temperature ($T_{\rm X})$ at small ($<R_{\rm e}$) and intermediate ($<5R_{\rm e}$) radii reasonably well. We find that the X-ray properties of lower mass galaxies depend on their star formation rates. In particular, in the magnitude range where the star-forming and quenched populations overlap, $M_{\rm K}\sim-24$ $ (M_*\sim10^{10.7}M_\odot)$, we find that the X-ray luminosities of star-forming galaxies are on average about an order of magnitude higher than those of their quenched counterparts. We show that this diversity in $L_{\rm X}$ is a direct manifestation of the quenching mechanism in the simulations, where the galaxies are quenched due to gas expulsion driven by SMBH kinetic feedback. The observed dichotomy in $L_{\rm X}$ is thus an important observable prediction for the SMBH feedback-based quenching mechanisms implemented in state-of-the-art cosmological simulations. While the current X-ray observations of star forming galaxies are broadly consistent with the predictions of the simulations, the observed samples are small and more decisive tests are expected from the sensitive all-sky X-ray survey with eROSITA.

We study the stellar-to-halo mass relation (SHMR) for central and satellite galaxies with total dynamical masses above 10^10.5 Msun using the suite of cosmological magneto-hydrodynamical simulations IllustrisTNG. In particular, we quantify environmental effects on satellite populations from TNG50, TNG100, and TNG300 located within the virial radius of group- and cluster-like hosts with total masses of 10^12-15.2 Msun. At fixed stellar mass, the satellite SHMR exhibits a distinct shift towards lower dynamical mass compared to the SHMR of centrals. Conversely, at fixed dynamical mass, satellite galaxies appear to have larger stellar-to-total mass fractions than centrals by up to a factor of a few. The systematic deviation from the central SHMR is larger for satellites in more massive hosts, at smaller cluster-centric distances, with earlier infall times, and that inhabit higher local density environments; moreover, it is in place already at early times (z < 2). Systematic environmental effects contribute to the perceived galaxy-to-galaxy variation in the measured SHMR when galaxies cannot be separated into satellites and centrals. The SHMR of satellites exhibits a larger scatter than centrals, over the whole range of dynamical mass (up to 0.8 dex). The shift of the satellite SHMR results mostly from tidal stripping of their dark matter, which affects satellites in an outside-in fashion: the departure of the satellite SHMR from the centrals' relation diminishes for measurements of dynamical mass in progressively smaller apertures. Finally, we provide a family of fitting functions for the SHMR predicted by IllustrisTNG.

Winds driven by stellar feedback are an essential part of the galactic ecosystem and are the main mechanism through which low-mass galaxies regulate their star formation. These winds are generally observed to be multi-phase with detections of entrained neutral and molecular gas. They are also thought to enrich the circum-galactic medium around galaxies with metals and dust. This ejected dust encodes information about the integrated star formation and outflow history of the galaxy. It is therefore, important to understand how much dust is entrained and driven out of the disc by galactic winds. Here we demonstrate that stellar feedback is efficient in driving dust-enriched winds and eject enough material to account for the amount of extraplanar dust observed in nearby galaxies. The amount of dust in the wind depends on the sites from where they are launched, with dustier galaxies launching more dust enriched outflows. Moreover, the outflowing cold-dense gas is significantly more dust-enriched than the volume filling hot tenuous material, naturally reproducing the complex multiphase structure of the outflowing wind observed in nearby galaxies. These results provide an important new insight into the dynamics, structure, and composition of galactic winds and their role in determining the dust content of the extragalactic gas in galaxies.

We investigate the fundamental plane (FP) evolution of early-type galaxies in the IllustrisTNG-100 simulation (TNG100) from redshift $z=0$ to $z=2$. We find that a tight plane relation already exists as early as $z=2$. Its scatter stays as low as $\sim 0.08$ dex across this redshift range. Both slope parameters $b$ and $c$ (where $R \propto \sigma^b I^c$ with $R$, $\sigma$, and $I$ being the typical size, velocity dispersion, and surface brightness) of the plane evolve mildly since $z=2$, roughly consistent with observations. The FP residual $\rm Res$ ($\equiv\,a\,+\,b\log \sigma\,+\,c\log I\,-\,\log R$, where $a$ is the zero point of the FP) is found to strongly correlate with stellar age, indicating that stellar age can be used as a crucial fourth parameter of the FP. However, we find that $4c+b+2=\delta$, where $\delta \sim 0.8$ for FPs in TNG, rather than zero as is typically inferred from observations. This implies that a tight power-law relation between the dynamical mass-to-light ratio $M_{\rm dyn}/L$ and the dynamical mass $M_{\rm dyn}$ (where $M_{\rm dyn}\equiv 5\sigma^2R/G$, with $G$ being the gravitational constant) is not present in the TNG100 simulation. Recovering such a relation requires proper mixing between dark matter and baryons, as well as star formation occurring with correct efficiencies at the right mass scales. This represents a powerful constraint on the numerical models, which has to be satisfied in future hydrodynamical simulations.

Surveys in the next decade will deliver large samples of galaxy clusters that transform our understanding of their formation. Cluster astrophysics and cosmology studies will become systematics limited with samples of this magnitude. With known properties, hydrodynamical simulations of clusters provide a vital resource for investigating potential systematics. However, this is only realized if we compare simulations to observations in the correct way. Here we introduce the \textsc{Mock-X} analysis framework, a multiwavelength tool that generates synthetic images from cosmological simulations and derives halo properties via observational methods. We detail our methods for generating optical, Compton-$y$ and X-ray images. Outlining our synthetic X-ray image analysis method, we demonstrate the capabilities of the framework by exploring hydrostatic mass bias for the IllustrisTNG, BAHAMAS and MACSIS simulations. Using simulation derived profiles we find an approximately constant bias $b\approx0.13$ with cluster mass, independent of hydrodynamical method or subgrid physics. However, the hydrostatic bias derived from synthetic observations is mass-dependent, increasing to $b=0.3$ for the most massive clusters. This result is driven by a single temperature fit to a spectrum produced by gas with a wide temperature distribution in quasi-pressure equilibrium. The spectroscopic temperature and mass estimate are biased low by cooler gas dominating the emission, due to its quadratic density dependence. The bias and the scatter in estimated mass remain independent of the numerical method and subgrid physics. Our results are consistent with current observations and future surveys will contain sufficient samples of massive clusters to confirm the mass dependence of the hydrostatic bias.

Recent cosmological hydrodynamical simulations are able to reproduce numerous statistical properties of galaxies that are consistent with observational data. Yet, the adopted subgrid models strongly affect the simulation outcomes, limiting the predictive power of these simulations. In this work, we perform a suite of isolated galactic disk simulations under the {\it SMUGGLE} framework and investigate how different subgrid models affect the properties of giant molecular clouds (GMCs). We employ {\sc astrodendro}, a hierarchical clump-finding algorithm, to identify GMCs in the simulations. We find that different choices of subgrid star formation efficiency, $\epsilon_{\rm ff}$, and stellar feedback channels, yield dramatically different mass and spatial distributions for the GMC populations. Without feedback, the mass function of GMCs has a shallower power-law slope and extends to higher mass ranges compared to runs with feedback. Moreover, higher $\epsilon_{\rm ff}$ results in faster molecular gas consumption and steeper mass function slopes. Feedback also suppresses power in the two-point correlation function (TPCF) of the spatial distribution of GMCs. Specifically, radiative feedback strongly reduces the TPCF on scales below 0.2~kpc, while supernova feedback reduces power on scales above 0.2~kpc. Finally, runs with higher $\epsilon_{\rm ff}$ exhibit a higher TPCF than runs with lower $\epsilon_{\rm ff}$, because the dense gas is depleted more efficiently thereby facilitating the formation of well-structured supernova bubbles. We argue that comparing simulated and observed GMC populations can help better constrain subgrid models in the next-generation of galaxy formation simulations.

We study the relation between the UV-slope, $\beta$, and the ratio between the IR and UV luminosities of galaxies, IRX, using TNG50, the latest installment of the IllustrisTNG galaxy formation simulations. We select 7280 TNG50 star-forming galaxies with stellar mass above $10^9M_\odot$ at selected redshifts, $0 \leq z \leq 4$, and perform radiative transfer calculations with SKIRT to model the effects of ISM dust on the simulated stellar light. We adopt a Milky Way (MW) type dust and a fixed dust-to-metal ratio of 0.3 throughout and find that TNG50 star-forming main-sequence galaxies agree with the empirically-derived reference IRX-$\beta$ relations at $z \lesssim 1$. There appears to be a redshift-dependent systematic offset with respect to the reference relations, with the TNG50 IRX-$\beta$ median relation steepening towards higher redshifts. This is to first order driven by variations in intrinsic UV-slopes due to different star-formation histories of galaxies selected at different cosmic epochs. Once differences in the intrinsic UV-slope are corrected for, TNG50 galaxies exhibit a scatter of 0.3 dex in IRX at fixed $\beta$ at all studied redshifts. These galaxy-to-galaxy variations are correlated with intrinsic galaxy properties, such that at fixed UV-slope galaxies with higher star--formation rates, star--formation efficiencies, gas metallicities, and stellar masses exhibit larger IRX values. We additionally demonstrate a degeneracy between stellar age and dust type. The combination of young stellar population age with a SMC type of dust for high-redshift ($z=4$) TNG50 galaxies can result in an IRX-$\beta$ relation similar to what is found for low-redshift galaxies with MW dust. This hampers the use of the IRX-$\beta$ relation as a proxy for dust type. We provide a redshift-dependent fitting function for the IRX-$\beta$ relation with MW dust based on our numerical models.

We explore the isothermal total density profiles of early-type galaxies (ETGs) in the IllustrisTNG simulation. For the selected 559 ETGs at $z = 0$ with stellar mass $10^{10.7}\mathrm{M}_{\odot} \leqslant M_{\ast} \leqslant 10^{11.9}\mathrm{M}_{\odot}$, the total power-law slope has a mean of $\langle\gamma^{\prime}\rangle = 2.011 \pm 0.007$ and a scatter of $\sigma_{\gamma^{\prime}} = 0.171$ over the radial range 0.4 to 4 times the stellar half mass radius. Several correlations between $\gamma^{\prime}$ and galactic properties including stellar mass, effective radius, stellar surface density, central velocity dispersion, central dark matter fraction and in-situ-formed stellar mass ratio are compared to observations and other simulations, revealing that IllustrisTNG reproduces many correlation trends, and in particular, $\gamma^{\prime}$ is almost constant with redshift below $z = 2$. Through analyzing IllustrisTNG model variations we show that black hole kinetic winds are crucial to lowering $\gamma^{\prime}$ and matching observed galaxy correlations. The effects of stellar winds on $\gamma^{\prime}$ are subdominant compared to AGN feedback, and differ due to the presence of AGN feedback from previous works. The density profiles of the ETG dark matter halos are well-described by steeper-than-NFW profiles, and they are steeper in the full physics (FP) run than their counterparts in the dark matter only (DMO) run. Their inner density slopes anti-correlates (remain constant) with the halo mass in the FP (DMO) run, and anti-correlates with the halo concentration parameter $c_{200}$ in both types of runs. The dark matter halos of low-mass ETGs are contracted whereas high-mass ETGs are expanded, suggesting that variations in the total density profile occur through the different halo responses to baryons.

Over the last decades, cosmological simulations of galaxy formation have been instrumental for advancing our understanding of structure and galaxy formation in the Universe. These simulations follow the non-linear evolution of galaxies modeling a variety of physical processes over an enormous range of scales. A better understanding of the physics relevant for shaping galaxies, improved numerical methods, and increased computing power have led to simulations that can reproduce a large number of observed galaxy properties. Modern simulations model dark matter, dark energy, and ordinary matter in an expanding space-time starting from well-defined initial conditions. The modeling of ordinary matter is most challenging due to the large array of physical processes affecting this matter component. Cosmological simulations have also proven useful to study alternative cosmological models and their impact on the galaxy population. This review presents a concise overview of the methodology of cosmological simulations of galaxy formation and their different applications.

We present a sample of 446 galaxy pairs constructed using the cosmological simulation IllustrisTNG-100 at z = 0, with M$_{FoF, dm}$ = 10$^{11}$-10$^{13.5}$ M$_{\odot}$. We produce ideal mock SDSS g-band images of all pairs to test the reliability of visual classification schema employed to produce samples of interacting galaxies. We visually classify each image as interacting or not based on the presence of a close neighbour, the presence of stellar debris fields, disturbed discs, and/or tidal features. By inspecting the trajectories of the pairs, we determine that these indicators correctly identify interacting galaxies $\sim45\%$ of the time. We subsequently split the sample into the visually identified interacting pairs (VIP; 38 pairs) and those which are interacting but are not visually identified (nonVIP; 47 pairs). We find that VIP have undergone a close passage nearly twice as recently as the nonVIP, and typically have higher stellar masses. Further, the VIP sit in dark matter haloes that are approximately 2.5 times as massive, in environments nearly 2 times as dense, and are almost a factor of 10 more affected by the tidal forces of their surroundings than the nonVIP. These factors conspire to increase the observability of tidal features and disturbed morphologies, making the VIP more likely to be identified. Thus, merger rate calculations which rely on stellar morphologies are likely to be significantly biased toward massive galaxy pairs which have recently undergone a close passage.

We present a model for the interaction between dust and radiation fields in the radiation hydrodynamic code AREPO-RT, which solves the moment-based radiative transfer equations on an unstructured moving mesh. Dust is directly treated using live simulation particles, each of which represent a population of grains that are coupled to hydrodynamic motion through a drag force. We introduce methods to calculate radiation pressure on and photon absorption by dust grains. By including a direct treatment of dust, we are able to calculate dust opacities and update radiation fields self-consistently based on the local dust distribution. This hybrid scheme coupling dust particles to an unstructured mesh for radiation is validated using several test problems with known analytic solutions, including dust driven via spherically-symmetric flux from a constant luminosity source and photon absorption from radiation incident on a thin layer of dust. Our methods are compatible with the multifrequency scheme in AREPO-RT, which treats UV and optical photons as single-scattered and IR photons as multi-scattered. At IR wavelengths, we model heating of and thermal emission from dust. Dust and gas are not assumed to be in local thermodynamic equilibrium but transfer energy through collisional exchange. We estimate dust temperatures by balancing these dust-radiation and dust-gas energy exchange rates. This framework for coupling dust and radiation can be applied in future radiation hydrodynamic simulations of galaxy formation.

The connections among galaxies, the dark matter halos where they form and the properties of the large-scale Cosmic Web still need to be completely disentangled. We use the cosmological hydrodynamical simulation TNG100 of the IllustrisTNG suite to quantify the effects played by the large-scale density field and the Cosmic Web morphology on the relation between halo mass and galaxy stellar mass. We select objects with total dynamical mass in the range $\geq 6.3\times 10^{10} h ^{-1}\, M_{\odot}$ up to a few $10^{14} h^{-1} \, M_{\odot}$ between redshift $z=4$ and redshift $z=0$. A Cosmic Web class (knot, filament, sheet, void) is assigned to each region of the volume using a density field deformation tensor-based method. We find that galaxy stellar mass strongly correlates with total dynamical mass and formation time, and more weakly with large-scale overdensity and Cosmic Web class. The latter two quantities correlate with each other, but are not entirely degenerate. Furthermore, we find that at fixed halo mass, galaxies with stellar mass lower than the median value are more likely to be found in voids and sheets, whereas galaxies with stellar mass higher than the median are more likely to be found in filaments and knots. Finally, we find that the dependence on environment is stronger for satellites than for centrals, and discuss the physical implications of these results.

We study the evolutionary trend of the total density profile of early-type galaxies (ETGs) in IllustrisTNG. To this end, we trace ETGs from $z=0$ to $z=4$ and measure the power-law slope $\gamma^{\prime}$ of the total density profile for their main progenitors. We find that their $\gamma^{\prime}$ steepen on average during $z\sim4-2$, then becoming shallower until $z=1$, after which they remain almost constant, aside from a residual trend of becoming shallower towards $z=0$. We also compare to a statistical sample of ETGs at different redshifts, selected based on their luminosity profiles and stellar masses. Due to different selection effects, the average slopes of the statistical samples follow a modified evolutionary trend. They monotonically decrease since $z=3$, and after $z\approx 1$, they remain nearly invariant with a mild increase towards $z=0$. These evolutionary trends are mass-dependent for both samples, with low-mass galaxies having in general steeper slopes than their more massive counterparts. Galaxies that transitioned to ETGs more recently have steeper mean slopes as they tend to be smaller and more compact at any given redshift. By analyzing the impact of mergers and AGN feedback on the progenitors' evolution, we conjecture a multi-phase path leading to isothermality in ETGs: dissipation associated with rapid wet mergers tends to steepen $\gamma^{\prime}$ from $z=4$ to $z=2$, whereas subsequent AGN feedback (especially in the kinetic mode) makes $\gamma^{\prime}$ shallower again from $z=2$ to $z=1$. Afterwards, passive evolution from $z=1$ to $z=0$, mainly through gas-poor mergers, mildly decreases $\gamma^{\prime}$ and maintains the overall mass distribution close to isothermal.

Bose-Einstein Condensate Dark Matter (BECDM; also known as Fuzzy Dark Matter) is motivated by fundamental physics and has recently received significant attention as a serious alternative to the established Cold Dark Matter (CDM) model. We perform cosmological simulations of BECDM gravitationally coupled to baryons and investigate structure formation at high redshifts ($z \gtrsim 5$) for a boson mass $m=2.5\cdot 10^{-22}~{\rm eV}$, exploring the dynamical effects of its wavelike nature on the cosmic web and the formation of first galaxies. Our BECDM simulations are directly compared to CDM as well as to simulations where the dynamical quantum potential is ignored and only the initial suppression of the power spectrum is considered -- a Warm Dark Matter-like ("WDM") model often used as a proxy for BECDM. Our simulations confirm that "WDM" is a good approximation to BECDM on large cosmological scales even in the presence of the baryonic feedback. Similarities also exist on small scales, with primordial star formation happening both in isolated haloes and continuously along cosmic filaments; the latter effect is not present in CDM. Global star formation and metal enrichment in these first galaxies are delayed in BECDM/"WDM" compared to the CDM case: in BECDM/"WDM" first stars form at $z\sim 13$/$13.5$ while in CDM star formation starts at $z\sim 35$. The signature of BECDM interference, not present in "WDM", is seen in the evolved dark matter power spectrum: although the small scale structure is initially suppressed, power on kpc scales is added at lower redshifts. Our simulations lay the groundwork for realistic simulations of galaxy formation in BECDM.

We present a novel framework to self-consistently model the effects of radiation fields, dust physics and molecular chemistry (H$_2$) in the interstellar medium (ISM) of galaxies. The model combines a state-of-the-art radiation hydrodynamics module with a non-equilibrium thermochemistry module that accounts for H$_2$ coupled to a realistic dust formation and destruction model, all integrated into the new stellar feedback framework SMUGGLE. We test this model on high-resolution isolated Milky-Way (MW) simulations. We show that photoheating from young stars makes stellar feedback more efficient, but this effect is quite modest in low gas surface density galaxies like the MW. The multi-phase structure of the ISM, however, is highly dependent on the strength of the interstellar radiation field. We are also able to predict the distribution of H$_2$, that allow us to match the molecular Kennicutt-Schmidt (KS) relation, without calibrating for it. We show that the dust distribution is a complex function of density, temperature and ionization state of the gas which cannot be reproduced by simple scaling relations often used in the literature. Our model is only able to match the observed dust temperature distribution if radiation from the old stellar population is considered, implying that these stars have a non-negligible contribution to dust heating in the ISM. Our state-of-the-art model is well-suited for performing next generation cosmological galaxy formation simulations, which will be able to predict a wide range of resolved ($\sim 10$ pc) properties of galaxies.

In hierarchical models of structure formation, the first galaxies form in low-mass dark matter potential wells, probing the behavior of dark matter on kiloparsec (kpc) scales. Even though these objects are below the detection threshold of current telescopes, future missions will open an observational window into this emergent world. In this Letter we investigate how the first galaxies are assembled in a `fuzzy' dark matter (FDM) cosmology where dark matter is an ultralight $\sim 10^{-22}$~eV boson and the primordial stars are expected to form along dense dark matter filaments. Using a first-of-its-kind cosmological hydrodynamical simulation, we explore the interplay between baryonic physics and unique wavelike features inherent to FDM. In our simulation, the dark matter filaments show coherent interference patterns on the boson de Broglie scale and develop cylindrical soliton-like cores which are unstable under gravity and collapse into kpc-scale spherical solitons. Features of the dark matter distribution are largely unaffected by the baryonic feedback. On the contrary, the distributions of gas and stars, which do form along the entire filament, exhibit central cores imprinted by dark matter -- a smoking gun signature of FDM.

We study black hole - host galaxy correlations, and the relation between the over-massiveness (the distance from the average $M_{BH}-\sigma$ relation) of super-massive black holes (SMBHs) and star formation histories of their host galaxies in the Illustris and TNG100 simulations. We find that both simulations are able to produce black hole scaling relations in general agreement with observations at $z=0$, but with noticeable discrepancies. Both simulations show an offset from the observations for the $M_{BH}-\sigma$ relation, and the relation between $M_{BH}$ and the Sersic index. The relation between $M_{BH}$ and stellar mass $M_*$ is tighter than the observations, especially for TNG100. For massive galaxies in both simulations, the hosts of over-massive SMBHs (those above the mean $M_{BH}-\sigma$ relation) tend to have larger Sersic indices and lower baryon conversion efficiency, suggesting a multidimensional link between SMBHs and properties of their hosts. In Illustris, the hosts of over-massive SMBHs have formed earlier and have lower present-day star formation rates, in qualitative agreement with the observations for massive galaxies with $\sigma>100 \rm km/s$. For low-mass galaxies, such a correlation still holds in Illustris but does not exist in the observed data. For TNG100, the correlation between SMBH over-massiveness and star formation history is much weaker. The hosts of over-massive SMBHs generally have consistently larger star formation rates throughout history. These galaxies have higher stellar mass as well, due to the strong $M_{BH}-M_*$ correlation. Our findings show that simulated SMBH scaling relations and correlations are sensitive to features in the modeling of SMBHs.

We present the star formation history (SFH) of the isolated (D~970 kpc) Local Group dwarf galaxy WLM measured from color-magnitude diagrams constructed from deep Hubble Space Telescope imaging. Our observations include a central (0.5 $r_h$) and outer field (0.7 $ r_h$) that reach below the oldest main sequence turnoff. WLM has no early dominant episode of star formation: 20% of its stellar mass formed by ~12.5 Gyr ago (z~5). It also has an SFR that rises to the present with 50% of the stellar mass within the most recent 5 Gyr (z<0.7). There is evidence of a strong age gradient: the mean age of the outer field is 5 Gyr older than the inner field despite being only 0.4 kpc apart. Some models suggest such steep gradients are associated with strong stellar feedback and dark matter core creation. The SFHs of real isolated dwarf galaxies and those from the the Feedback In Realistic Environment suite are in good agreement for $M_{\star}(z=0) \sim 10^7-10^9 M_{\odot}$, but in worse agreement at lower masses ($M_{\star}(z=0) \sim 10^5-10^7 M_{\odot}$). These differences may be explainable by systematics in the models (e.g., reionization model) and/or observations (HST field placement). We suggest that a coordinated effort to get deep CMDs between HST/JWST (crowded central fields) and WFIRST (wide-area halo coverage) is the optimal path for measuring global SFHs of isolated dwarf galaxies.

We present a new cosmological, magnetohydrodynamical simulation for galaxy formation: TNG50, the third and final installment of the IllustrisTNG project. TNG50 evolves 2x2160^3 dark-matter particles and gas cells in a volume 50 comoving Mpc across. It hence reaches a numerical resolution typical of zoom-in simulations, with a baryonic element mass of 8.5x10^4 Msun and an average cell size of 70-140 parsecs in the star-forming regions of galaxies. Simultaneously, TNG50 samples ~700 (6,500) galaxies with stellar masses above 10^10 (10^8) Msun at z=1. Here we investigate the structural and kinematical evolution of star-forming galaxies across cosmic time (0 < z < 6). We quantify their sizes, disk heights, 3D shapes, and degree of rotational vs. dispersion-supported motions as traced by rest-frame V-band light (i.e. roughly stellar mass) and by Halpha light (i.e. star-forming and dense gas). The unprecedented resolution of TNG50 enables us to model galaxies with sub-kpc half-light radii and with <300-pc disk heights. Coupled with the large-volume statistics, we characterize a diverse, redshift- and mass-dependent structural and kinematical morphological mix of galaxies all the way to early epochs. Our model predicts that for star-forming galaxies the fraction of disk-like morphologies, based on 3D stellar shapes, increases with both cosmic time and galaxy stellar mass. Gas kinematics reveal that the vast majority of 10^9-11.5 Msun star-forming galaxies are rotationally-supported disks for most cosmic epochs (Vmax/sigma>2-3, z<5), being dynamically hotter at earlier epochs (z>1.5). Despite large velocity dispersion at high redshift, cold and dense gas in galaxies predominantly arranges in disky or elongated shapes at all times and masses; these gaseous components exhibit rotationally-dominated motions far exceeding the collisionless stellar bodies.

We introduce the Stars and MUltiphase Gas in GaLaxiEs -- SMUGGLE model, an explicit and comprehensive stellar feedback model for the moving-mesh code arepo. This novel sub-resolution model resolves the multiphase gas structure of the interstellar medium and self-consistently generates gaseous outflows. The model implements crucial aspects of stellar feedback including photoionization, radiation pressure, energy and momentum injection from stellar winds and from supernovae. We explore this model in high-resolution isolated simulations of Milky Way-like disc galaxies. Stellar feedback regulates star formation to the observed level and naturally captures the establishment of a Kennicutt-Schmidt relation. This result is achieved independent of the numerical mass and spatial resolution of the simulations. Gaseous outflows are generated with average mass loading factors of the order of unity. Strong outflow activity is correlated with peaks in the star formation history of the galaxy with evidence that most of the ejected gas eventually rains down onto the disc in a galactic fountain flow that sustains late-time star formation. Finally, the interstellar gas in the galaxy shows a distinct multiphase distribution with a coexistence of cold, warm and hot phases.

Milky Way (MW) satellites reside within dark matter (DM) subhalos with a broad distribution of circular velocity profiles. This diversity is enhanced with the inclusion of ultra-faint satellites, which seemingly have very high DM densities, albeit with large systematic uncertainties. We argue that if confirmed, this large diversity in the MW satellite population poses a serious test for the structure formation theory with possible implications for the DM nature. For the Cold Dark Matter model, the diversity might be a signature of the combined effects of subhalo tidal disruption by the MW disk and strong supernova feedback. For models with a dwarf-scale cutoff in the power spectrum, the diversity is a consequence of the lower abundance of dwarf-scale halos. This diversity is most challenging for Self-Interacting Dark Matter (SIDM) models with cross sections $\sigma/m_\chi\gtrsim1~$cm$^2$g$^{-1}$ where subhalos have too low densities to explain the ultra-faint galaxies. We propose a novel solution to explain the diversity of MW satellites based on the gravothermal collapse of SIDM haloes. This solution requires a velocity-dependent cross section that predicts a bimodal distribution of cuspy dense (collapsed) subhaloes consistent with the ultra-faint satellites, and cored lower density subhaloes consistent with the brighter satellites.

We present the new TNG50 cosmological, magnetohydrodynamical simulation -- the third and final volume of the IllustrisTNG project. This simulation occupies a unique combination of large volume and high resolution, with a 50 Mpc box sampled by 2160^3 gas cells (baryon mass of 8x10^4 Msun). The median spatial resolution of star-forming ISM gas is ~100-140 parsecs. This resolution approaches or exceeds that of modern 'zoom' simulations of individual massive galaxies, while the volume contains ~20,000 resolved galaxies with M*>10^7 Msun. Herein we show first results from TNG50, focusing on galactic outflows driven by supernovae as well as supermassive black hole feedback. We find that the outflow mass loading is a non-monotonic function of galaxy stellar mass, turning over and rising rapidly above 10^10.5 Msun due to the action of the central black hole. Outflow velocity increases with stellar mass, and at fixed mass is faster at higher redshift. The TNG model can produce high velocity, multi-phase outflows which include cool, dense components. These outflows reach speeds in excess of 3000 km/s out to 20 kpc with an ejective, BH-driven origin. Critically, we show how the relative simplicity of model inputs (and scalings) at the injection scale produces complex behavior at galactic and halo scales. For example, despite isotropic wind launching, outflows exhibit natural collimation and an emergent bipolarity. Furthermore, galaxies above the star-forming main sequence drive faster outflows, although this correlation inverts at high mass with the onset of quenching, whereby low luminosity, slowly accreting, massive black holes drive the strongest outflows.

We provide an overview of the science case, instrument configuration and project plan for the next-generation ground-based cosmic microwave background experiment CMB-S4, for consideration by the 2020 Decadal Survey.

We present the science case, reference design, and project plan for the Stage-4 ground-based cosmic microwave background experiment CMB-S4.

Simulating the dust content of galaxies and their surrounding gas is challenging due to the wide range of physical processes affecting the dust evolution. Here we present cosmological hydrodynamical simulations of a cluster of galaxies, $M_\text{200,crit}=6 \times 10^{14}\,{\rm M_\odot}$, including a novel dust model for the moving mesh code {\sc Arepo}. This model includes dust production, growth, supernova-shock-driven destruction, ion-collision-driven thermal sputtering, and high temperature dust cooling through far infrared re-radiation of collisionally deposited electron energies. Adopting a rather low thermal sputtering rate, we find, consistent with observations, a present-day overall dust-to-gas ratio of $\sim 2\times 10^{-5}$, a total dust mass of $\sim 2\times 10^9\,{\rm M_\odot}$, and a dust mass fraction of $\sim 3\times 10^{-6}$. The typical thermal sputtering timescales within $\sim 100\,{\rm kpc}$ are around $\sim 10\,{\rm Myr}$, and increase towards the outer parts of the cluster to $\sim 10^3\,{\rm Myr}$ at a cluster-centric distance of $1\,{\rm Mpc}$. The condensation of gas phase metals into dust grains reduces high temperature metal-line cooling, but also leads to additional dust infrared cooling. The additional infrared cooling changes the overall cooling rate in the outer parts of the cluster, beyond $\sim 1\,{\rm Mpc}$, by factors of a few. This results in noticeable changes of the entropy, temperature, and density profiles of cluster gas once dust formation is included. The emitted dust infrared emission due to dust cooling is consistent with observational constraints.

Using the IllustrisTNG simulations, we investigate the connection between galaxy morphology and star formation in central galaxies with stellar masses in the range $10^9-10^{11.5}~\mathrm{M}_{\odot}$. We quantify galaxy morphology by a kinematical decomposition of the stellar component into a spheroidal and a disc component (spheroid-to-total ratio, S/T) and by the concentration of the stellar mass density profile ($C_{82}$). S/T is correlated with stellar mass and star-formation activity, while $C_{82}$ correlates only with stellar mass. Overall, we find good agreement with observational estimates for both S/T and $C_{82}$. Low and high mass galaxies are dominated by random stellar motion, while only intermediate-mass galaxies ($M_{\star}\approx10^{10}-10^{10.5}~\mathrm{M}_{\odot}$) are dominated by ordered rotation. Whereas higher-mass galaxies are typical spheroids with high concentrations, lower-mass galaxies have low concentration, pointing to different formation channels. Although we find a correlation between S/T and star-formation activity, in the TNG model galaxies do not necessarily change their morphology when they transition through the green valley or when they cease their star formation, this depending on galaxy stellar mass and morphological estimator. Instead, the morphology (S/T and $C_{82}$) is generally set during the star-forming phase of galaxies. The apparent correlation between S/T and star formation arises because earlier-forming galaxies had, on average, a higher S/T at a given stellar mass. Furthermore, we show that mergers drive in-situ bulge formation in intermediate-mass galaxies and are responsible for the recent spheroidal mass assembly in the massive galaxies with $M_{\star}>10^{11}~\mathrm{M}_{\odot}$. In particular, these massive galaxies assemble about half of the spheroidal mass while star-forming and the other half through mergers while quiescent.

We present self-consistent radiation hydrodynamic simulations of hydrogen reionization performed with Arepo-RT complemented by a state-of-the-art galaxy formation model. We examine how photoheating feedback, due to reionization, shapes the galaxies properties. Our fiducial model completes reionization by $z\approx6$ and matches observations of the Ly$\alpha$ forest, the CMB electron scattering optical depth, the high-redshift UV luminosity function, and stellar mass function. Contrary to previous works, photoheating suppresses star formation rates by more than $50\%$ only in halos less massive than $\sim10^{8.4}\ M_\odot$ ($\sim10^{8.8}\ M_\odot$) at $z=6$ $(z=5)$, suggesting inefficient photoheating feedback from photons within galaxies. The use of a uniform UV background that heats up the gas at $z\approx10.7$ generates an earlier onset of suppression of star formation compared to our fiducial model. This discrepancy can be mitigated by adopting a UV background model with a more realistic reionization history. In the absence of stellar feedback, photoheating alone is only able to quench halos less massive than $\sim10^9\ M_\odot$ at $z\gtrsim5$, implying that photoheating feedback is sub-dominant in regulating star formation. In addition, stellar feedback, implemented as a non-local galactic wind scheme in the simulations, weakens the strength of photoheating feedback by reducing the amount of stellar sources. Most importantly, photoheating does not leave observable imprints in the UV luminosity function, stellar mass function, or the cosmic star formation rate density. The feasibility of using these observables to detect imprints of reionization therefore requires further investigation.

We have recently developed a post-processing framework to estimate the abundance of atomic and molecular hydrogen (HI and H2, respectively) in galaxies in large-volume cosmological simulations. Here we compare the HI and H2 content of IllustrisTNG galaxies to observations. We mostly restrict this comparison to $z \approx 0$ and consider six observational metrics: the overall abundance of HI and H2, their mass functions, gas fractions as a function of stellar mass, the correlation between H2 and star formation rate, the spatial distribution of gas, and the correlation between gas content and morphology. We find generally good agreement between simulations and observations, particularly for the gas fractions and the HI mass-size relation. The H2 mass correlates with star formation rate as expected, revealing an almost constant depletion time that evolves up to z = 2 as observed. However, we also discover a number of tensions with varying degrees of significance, including an overestimate of the total neutral gas abundance at z = 0 by about a factor of two and a possible excess of satellites with no or very little neutral gas. These conclusions are robust to the modelling of the HI/H2 transition. In terms of their neutral gas properties, the IllustrisTNG simulations represent an enormous improvement over the original Illustris run. All data used in this paper are publicly available as part of the IllustrisTNG data release.

We present a machine-learning approach for estimating galaxy cluster masses from Chandra mock images. We utilize a Convolutional Neural Network (CNN), a deep machine learning tool commonly used in image recognition tasks. The CNN is trained and tested on our sample of 7,896 Chandra X-ray mock observations, which are based on 329 massive clusters from the IllustrisTNG simulation. Our CNN learns from a low resolution spatial distribution of photon counts and does not use spectral information. Despite our simplifying assumption to neglect spectral information, the resulting mass values estimated by the CNN exhibit small bias in comparison to the true masses of the simulated clusters (-0.02 dex) and reproduce the cluster masses with low intrinsic scatter, 8% in our best fold and 12% averaging over all. In contrast, a more standard core-excised luminosity method achieves 15-18% scatter. We interpret the results with an approach inspired by Google DeepDream and find that the CNN ignores the central regions of clusters, which are known to have high scatter with mass.

Supermassive black hole feedback is thought to be responsible for the lack of star formation, or quiescence, in a significant fraction of galaxies. We explore how observable correlations between the specific star formation rate (sSFR), stellar mass (M$_{\rm{star}}$), and black hole mass (M$_{\rm{BH}}$) are sensitive to the physics of black hole feedback in a galaxy formation model. We use the IllustrisTNG simulation suite, specifically the TNG100 simulation and ten model variations that alter the parameters of the black hole model. Focusing on central galaxies at $z = 0$ with M$_{\rm{star}} > 10^{10}$ M$_{\odot}$, we find that the sSFR of galaxies in IllustrisTNG decreases once the energy from black hole kinetic winds at low accretion rates becomes larger than the gravitational binding energy of gas within the galaxy stellar radius. This occurs at a particular M$_{\rm{BH}}$ threshold above which galaxies are found to sharply transition from being mostly star-forming to mostly quiescent. As a result of this behavior, the fraction of quiescent galaxies as a function of M$_{\rm{star}}$ is sensitive to both the normalization of the M$_{\rm{BH}}$-M$_{\rm{star}}$ relation and the M$_{\rm{BH}}$ threshold for quiescence in IllustrisTNG. Finally, we compare these model results to observations of 91 central galaxies with dynamical M$_{\rm{BH}}$ measurements with the caveat that this sample is not representative of the whole galaxy population. While IllustrisTNG reproduces the observed trend that quiescent galaxies host more massive black holes, the observations exhibit a broader scatter in M$_{\rm{BH}}$ at a given M$_{\rm{star}}$ and show a smoother decline in sSFR with M$_{\rm{BH}}$.

We perform a series of cosmological hydrodynamic simulations to investigate the effects of non-gravitational dark matter (DM) interactions on the intergalactic medium (IGM). In particular, we use the ETHOS framework (Cyr-Racine et al. 2016; Vogelsberger et al. 2016) to compare statistics of the Lyman-$\alpha$ forest in cold dark matter (CDM) with an alternative model in which the DM couples strongly with a relativistic species in the early universe. These models are characterised by a cutoff in the linear power spectrum, followed by a series of 'dark acoustic oscillations' (DAOs) on sub-dwarf scales. While the primordial cutoff delays the formation of the first galaxies, structure builds-up more rapidly in the interacting DM model compared to CDM. We show that although DAOs are quickly washed away in the non-linear clustering of DM at $z\lesssim10$, their signature can be imprinted prominently in the Lyman-$\alpha$ flux power spectrum at $z>5$. On scales larger than the cutoff ($k\sim0.08$ s/km for the specific model considered here), the relative difference to CDM is reminiscent of a warm dark matter (WDM) model with a similar initial cutoff; however, the redshift evolution on smaller scales is distinctly different. The appearance and disappearance of DAOs in the Lyman-$\alpha$ flux spectrum provides a powerful way to distinguish interacting DM models from WDM and, indeed, variations in the thermal history of the IGM that may also induce a small-scale cutoff.

Motivated by the recently discovered kinematic "Hubble sequence" shown by the stellar orbit-circularity distribution of 260 CALIFA galaxies, we make use of a comparable galaxy sample at z = 0 with a stellar mass range from 5E9 to 5E11 solar masses, selected from the IllustrisTNG simulation and study their stellar orbit compositions in relation to a number of other fundamental galaxy properties.We find that the TNG100 simulation broadly reproduces the observed fractions of different orbital components and their stellar mass dependencies. In particular, the mean mass dependencies of the luminosity fractions for the kinematically warm and hot orbits are well reproduced within model uncertainties of the observed galaxies. The simulation also largely reproduces the observed peak and trough features at a stellar mass of 1-2E10 solar masses, in the mean distributions of the cold- and hot-orbit fractions, respectively, indicating fewer cooler orbits and more hotter orbits in both more- and less-massive galaxies beyond such a mass range. Several marginal disagreements are seen between the simulation and observations: the average cold-orbit (counter-rotating) fractions of the simulated galaxies below (above) a stellar mass of 6E10 solar masses, are systematically higher than the observational data by < 10% (absolute orbital fraction); the simulation also seems to produce more scatter for the cold-orbit fraction and less so for the non-cold orbits at any given galaxy mass. Possible causes that stem from the adopted heating mechanisms are discussed.

Supersonically induced gas objects (SIGOs) with little to no dark matter component are predicted to exist in patches of the Universe with non-negligible relative velocity between baryons and the dark matter at the time of recombination. Using {\sc arepo} hydrodynamic simulations we find that the gas densities inside these objects are high enough to allow stars to form. An estimate of the luminosity of the first star clusters formed within these SIGOs suggests that they may be observed at high redshift using future HST and JWST observations. Furthermore, our simulations indicate that SIGOs lie in a distinct place in the luminosity-radius parameter space, which can be used observationally to distinguish SIGOs from dark-matter hosting gas systems. Finally, as a proof-of-concept, we model star formation before reionization and evolve these systems to current times. We find that SIGOs occupy a similar part of the magnitude-radius parameter space as globular clusters. These results suggest that SIGOs may be linked with present-day metal-poor local globular clusters. Since the relative velocity between the baryons and dark matter is coherent over a few Mpc scales, we predict that if this is the dominant mechanism for the formation of globular clusters, their abundance should vary significantly over these scales.

We use the IllustrisTNG (TNG) simulations to explore the galaxy-halo connection as inferred from state-of-the-art cosmological, magnetohydrodynamical simulations. With the high mass resolution and large volume achieved by combining the 100 Mpc (TNG100) and 300 Mpc (TNG300) volumes, we establish the mean occupancy of central and satellite galaxies and their dependence on the properties of the dark matter haloes hosting them. We derive best-fitting HOD parameters from TNG100 and TNG300 for target galaxy number densities of $\bar{n}_g = 0.032\,h^3$Mpc$^{-3}$ and $\bar{n}_g = 0.016\,h^3$Mpc$^{-3}$, respectively, corresponding to a minimum galaxy stellar mass of $M_\star\sim1.9\times10^9\,{\rm M}_\odot$ and $M_\star\sim3.5\times10^9\,{\rm M}_\odot$, respectively, in hosts more massive than $10^{11}\,{\rm M}_\odot$. Consistent with previous work, we find that haloes located in dense environments, with low concentrations, later formation times, and high angular momenta are richest in their satellite population. At low mass, highly-concentrated haloes and those located in overdense regions are more likely to contain a central galaxy. The degree of environmental dependence is sensitive to the definition adopted for the physical boundary of the host halo. We examine the extent to which correlations between galaxy occupancy and halo properties are independent and demonstrate that HODs predicted by halo mass and present-day concentration capture the qualitative dependence on the remaining halo properties. At fixed halo mass, concentration is a strong predictor of the stellar mass of the central galaxy, which may play a defining role in the fate of the satellite population. The radial distribution of satellite galaxies, which exhibits a universal form across a wide range of host halo mass, is described accurately by the best-fit NFW density profile of their host haloes.

Energetic feedback from star clusters plays a pivotal role in shaping the dynamical evolution of giant molecular clouds (GMCs). To study the effects of stellar feedback on the star formation efficiency of the clouds and the dynamical response of embedded star clusters, we perform a suite of isolated GMC simulations with star formation and momentum feedback subgrid models using the moving-mesh hydrodynamics code \textsc{Arepo}. The properties of our simulated GMCs span a wide range of initial mass, radius, and velocity configurations. We find that the ratio of the final stellar mass to the total cloud mass, $\epsilon_{\rm int}$, scales strongly with the initial cloud surface density and momentum feedback strength. This correlation is explained by an analytic model that considers force balancing between gravity and momentum feedback. For all simulated GMCs, the stellar density profiles are systematically steeper than that of the gas at the epochs of the peaks of star formation, suggesting a centrally concentrated stellar distribution. We also find that star clusters are always in a sub-virial state with a virial parameter $\sim0.6$ prior to gas expulsion. Both the sub-virial dynamical state and steeper stellar density profiles prevent clusters from dispersal during the gas removal phase of their evolution. The final cluster bound fraction is a continuously increasing function of $\epsilon_{\rm int}$. GMCs with star formation efficiency smaller than 0.5 are still able to form clusters with large bound fractions.

We analyze the IllustrisTNG simulations to study the mass, volume fraction and phase distribution of gaseous baryons embedded in the knots, filaments, sheets and voids of the Cosmic Web from redshift $z=8$ to redshift $z=0$. We find that filaments host more star-forming gas than knots, and that filaments also have a higher relative mass fraction of gas in this phase than knots. We also show that the cool, diffuse Intergalactic Medium (IGM; $T<10^5 \, {\rm K}$, $ n_{\rm H}<10^{-4}(1+z) \, {\rm cm^{-3}}$) and the Warm-Hot Intergalactic Medium (WHIM; $ 10^5 \, {\rm K} <T<10^7 \, {\rm K}$, $ n_{\rm H} <10^{-4}(1+z)\, {\rm cm^{-3}}$) constitute $\sim 39\%$ and $\sim 46\%$ of the baryons at redshift $z=0$, respectively. Our results indicate that the WHIM may constitute the largest reservoir of {\it missing} baryons at redshift $z=0$. Using our Cosmic Web classification, we predict the WHIM to be the dominant baryon mass contribution in filaments and knots at redshift $z=0$, but not in sheets and voids where the cool, diffuse IGM dominates. We also characterise the evolution of WHIM and IGM from redshift $z=4$ to redshift $z=0$, and find that the mass fraction of WHIM in filaments and knots evolves only by a factor $\sim 2$ from redshift $z=0$ to $z=1$, but declines faster at higher redshift. The WHIM only occupies $4-11\%$ of the volume at redshift $0\leq z \leq 1$. We predict the existence of a significant number of currently undetected OVII and NeIX absorption systems in cosmic filaments which could be detected by future X-ray telescopes like Athena.

The James Webb Space Telescop (JWST) promises to revolutionise our understanding of the early Universe, and contrasting its upcoming observations with predictions of the $\Lambda$CDM model requires detailed theoretical forecasts. Here, we exploit the large dynamic range of the IllustrisTNG simulation suite, TNG50, TNG100, and TNG300, to derive multi-band galaxy luminosity functions from $z=2$ to $z=10$. We put particular emphasis on the exploration of different dust attenuation models to determine galaxy luminosity functions for the rest-frame ultraviolet (UV), and apparent wide NIRCam bands. Our most detailed dust model is based on continuum Monte Carlo radiative transfer calculations employing observationally calibrated dust properties. This calibration results in constraints on the redshift evolution of the dust attenuation normalisation and dust-to-metal ratios yielding a stronger redshift evolution of the attenuation normalisation compared to most previous theoretical studies. Overall we find good agreement between the rest-frame UV luminosity functions and observational data for all redshifts, also beyond the regimes used for the dust-model calibrations. Furthermore, we also recover the observed high redshift ($z=4-6$) UV luminosity versus stellar mass relation, the H$\alpha$ versus star formation rate relation, and the H$\alpha$ luminosity function at $z=2$. The bright end ($M_{\rm UV}>-19.5$) cumulative galaxy number densities are consistent with observational data. For the F200W NIRCam band, we predict that JWST will detect $\sim 80$ ($\sim 200$) galaxies with a signal-to-noise ratio of $10$ ($\sim 5$) within the NIRCam field of view, $2.2\times2.2 \,{\rm arcmin}^{2}$, for a total exposure time of $10^5{\rm s}$ in the redshift range $z=8 \pm 0.5$. These numbers drop to $\sim 10$ ($\sim 40$) for an exposure time of $10^4{\rm s}$.

We present image-based evolution of galaxy mergers from the Illustris cosmological simulation at 12 time-steps over 0.5 < z < 5. To do so, we created approximately one million synthetic deep Hubble Space Telescope and James Webb Space Telescope images and measured common morphological indicators. Using the merger tree, we assess methods to observationally select mergers with stellar mass ratios as low as 10:1 completing within +/- 250 Myr of the mock observation. We confirm that common one- or two-dimensional statistics select mergers so defined with low purity and completeness, leading to high statistical errors. As an alternative, we train redshift-dependent random forests (RFs) based on 5-10 inputs. Cross-validation shows the RFs yield superior, yet still imperfect, measurements of the late-stage merger fraction, and they select more mergers in bulge-dominated galaxies. When applied to CANDELS morphology catalogs, the RFs estimate a merger rate increasing to at least z = 3, albeit two times higher than expected by theory. This suggests possible mismatches in the feedback-determined morphologies, but affirms the basic understanding of galaxy merger evolution. The RFs achieve completeness of roughly 70% at 0.5 < z < 3, and purity increasing from 10% at z = 0.5 to 60% at z = 3. At earlier times, the training sets are insufficient, motivating larger simulations and smaller time sampling. By blending large surveys and large simulations, such machine learning techniques offer a promising opportunity to teach us the strengths and weaknesses of inferences about galaxy evolution.

We present the full public release of all data from the TNG100 and TNG300 simulations of the IllustrisTNG project. IllustrisTNG is a suite of large volume, cosmological, gravo-magnetohydrodynamical simulations run with the moving-mesh code Arepo. TNG includes a comprehensive model for galaxy formation physics, and each TNG simulation self-consistently solves for the coupled evolution of dark matter, cosmic gas, luminous stars, and supermassive blackholes from early time to the present day, z=0. Each of the flagship runs -- TNG50, TNG100, and TNG300 -- are accompanied by lower-resolution and dark-matter only counterparts, and we discuss scientific and numerical cautions and caveats relevant when using TNG. Full volume snapshots are available at 100 redshifts; halo and subhalo catalogs at each snapshot and merger trees are also released. The data volume now directly accessible online is ~750 TB, including 1200 full volume snapshots and ~80,000 high time-resolution subbox snapshots. This will increase to ~1.1 PB with the future release of TNG50. Data access and analysis examples are available in IDL, Python, and Matlab. We describe improvements and new functionality in the web-based API, including on-demand visualization and analysis of galaxies and halos, exploratory plotting of scaling relations and other relationships between galactic and halo properties, and a new JupyterLab interface. This provides an online, browser-based, near-native data analysis platform which supports user computation with fully local access to TNG data, alleviating the need to download large simulated datasets.

The ALMA Spectroscopic Survey in the Hubble Ultra Deep Field (ASPECS) provides new constraints for galaxy formation models on the molecular gas properties of galaxies. We compare results from ASPECS to predictions from two cosmological galaxy formation models: the IllustrisTNG hydrodynamical simulations and the Santa Cruz semi-analytic model (SC SAM). We explore several recipes to model the H$_2$ content of galaxies, finding them to be consistent with one another, and take into account the sensitivity limits and survey area of ASPECS. For a canonical CO-to-H$_2$ conversion factor of $\alpha_{\rm CO} = 3.6\,\rm{M}_\odot/(\rm{K}\,\rm{km/s}\,\rm{pc}^{2})$ the results of our work include: (1) the H$_2$ mass of $z>1$ galaxies predicted by the models as a function of their stellar mass is a factor of 2-3 lower than observed; (2) the models do not reproduce the number of H$_2$-rich ($M_{\rm H2} > 3\times 10^{10}\,\rm{M}_\odot$) galaxies observed by ASPECS; (3) the H$_2$ cosmic density evolution predicted by IllustrisTNG (the SC SAM) is in tension (only just agrees) with the observed cosmic density, even after accounting for the ASPECS selection function and field-to-field variance effects. The tension between models and observations at $z>1$ can be alleviated by adopting a CO-to-H$_2$ conversion factor in the range $\alpha_{\rm CO} = 2.0 - 0.8\,\rm{M}_\odot/(\rm{K}\,\rm{km/s}\,\rm{pc}^{2})$. Additional work on constraining the CO-to-H$_2$ conversion factor and CO excitation conditions of galaxies through observations and theory will be necessary to more robustly test the success of galaxy formation models.

We analyze the optical morphologies of galaxies in the IllustrisTNG simulation at $z\sim0$ with a Convolutional Neural Network trained on visual morphologies in the Sloan Digital Sky Survey. We generate mock SDSS images of a mass complete sample of $\sim12,000$ galaxies in the simulation using the radiative transfer code SKIRT and include PSF and noise to match the SDSS r-band properties. The images are then processed through the exact same neural network used to estimate SDSS morphologies to classify simulated galaxies in four morphological classes (E, S0/a, Sab, Scd). The CNN model finds that $\sim95\%$ of the simulated galaxies fall in one the four main classes with high confidence. The mass-size relations of the simulated galaxies divided by morphological type also reproduce well the slope and the normalization of observed relations which confirms the realism of optical morphologies in the TNG suite. However, the Stellar Mass Functions decomposed into different morphologies still show significant discrepancies with observations both at the low and high mass end. We find that the high mass end of the SMF is dominated in TNG by massive disk galaxies while early-type galaxies dominate in the observations according to the CNN classifications. The present work highlights the importance of detailed comparisons between observations and simulations in comparable conditions.

Self-interacting dark matter provides a promising alternative for the cold dark matter paradigm to solve potential small-scale galaxy formation problems. Nearly all self-interacting dark matter simulations so far have considered only elastic collisions. Here we present simulations of a galactic halo within a generic inelastic model using a novel numerical implementation in the Arepo code to study arbitrary multi-state inelastic dark matter scenarios. For this model we find that inelastic self-interactions can: (i) create larger subhalo density cores compared to elastic models for the same cross section normalisation; (ii) lower the abundance of satellites without the need for a power spectrum cutoff; (iii) reduce the total halo mass by about 10%; (iv) inject the energy equivalent of O(100) million Type II supernovae in galactic haloes through level de-excitation; (v) avoid the gravothermal catastrophe due to removal of particles from halo centers. We conclude that a ~5 times larger elastic cross section is required to achieve the same central density reduction as the inelastic model. This implies that well-established constraints on self-interacting cross sections have to be revised if inelastic collisions are the dominant mode. In this case significantly smaller cross sections can achieve the same core density reduction thereby increasing the parameter space of allowed models considerably.

We select galaxies from the IllustrisTNG hydrodynamical simulations ($M_*>10^9~\rm M_\odot$ at $0\le z\le2$) and characterize the shapes and evolutions of their UVJ and star-formation rate -- stellar mass (SFR-$M_*$) diagrams. We quantify the systematic uncertainties related to different criteria to classify star-forming vs. quiescent galaxies, different SFR estimates, and by accounting for the star formation measured within different physical apertures. The TNG model returns the observed features of the UVJ diagram at $z\leq2$, with a clear separation between two classes of galaxies. It also returns a tight star-forming main sequence (MS) for $M_*<10^{10.5}\,\rm M_\odot$ with a $\sim0.3$ dex scatter at $z\sim0$ in our fiducial choices. If a UVJ-based cut is adopted, the TNG MS exhibits a downwardly bending at stellar masses of about $10^{10.5-10.7}~\rm M_\odot$. Moreover, the model predicts that $\sim80\,(50)$ per cent of $10^{10.5-11}~\rm M_\odot$ galaxies at $z=0~(z=2)$ are quiescent and the numbers of quenched galaxies at intermediate redshifts and high masses are in better agreement with observational estimates than previous models. However, shorter SFR-averaging timescales imply higher normalizations and scatter of the MS, while smaller apertures lead to underestimating the galaxy SFRs: overall we estimate the inspected systematic uncertainties to sum up to about $0.2-0.3$ dex in the locus of the MS and to about 15 percentage points in the quenched fractions. While TNG color distributions are clearly bimodal, this is not the case for the SFR logarithmic distributions in bins of stellar mass (SFR$\geq 10^{-3}~\rm M_\odot$yr$^{-1}$). Finally, the slope and $z=0$ normalization of the TNG MS are consistent with observational findings; however, the locus of the TNG MS remains lower by about $0.2-0.5$ dex at $0.75\le z<2$ than the available observational estimates taken at face value.

We study the abundance of subhaloes in the hydrodynamical cosmological simulation Illustris, which includes both baryons and dark matter in a $\Lambda$CDM volume 106.5 Mpc a side. We compare Illustris to its dark matter-only (DMO) analogue, Illustris-Dark, and quantify the effects of baryonic processes on the demographics of subhaloes in the host mass range $10^{11}$ to $3 \times 10^{14} \msun$. We focus on both the evolved ($z=0$) subhalo cumulative mass functions (SHMF) and the statistics of subhaloes ever accreted, i.e. infall subhalo mass function. We quantify the variance in subhalo abundance at fixed host mass and investigate the physical reasons responsible for such scatter. We find that in Illustris, baryonic physics impacts both the infall and $z=0$ subhalo abundance by tilting the DMO function and suppressing the abundance of low-mass subhaloes. The breaking of self-similarity in the subhalo abundance at $z=0$ is enhanced by the inclusion of baryonic physics. The non-monotonic alteration of the evolved subhalo abundances can be explained by the modification of the concentration--mass relation of Illustris hosts compared to Illustris-Dark.Interestingly, the baryonic implementation in Illustris does not lead to an increase in the halo-to-halo variation compared to Illustris-Dark. In both cases, the normalized intrinsic scatter today is larger for Milky Way-like haloes than for cluster-sized objects. For Milky Way-like haloes, it increases from about eight per cent at infall to about 25 per cent at the current epoch. In both runs, haloes of fixed mass formed later host more subhaloes than early formers.

Achieving a precise understanding of galaxy formation in a cosmological context is one of the great challenges in theoretical astrophysics, due to the vast range of spatial scales involved in the relevant physical processes. Observations in the millimeter bands, particularly those using the cosmic microwave background (CMB) radiation as a "backlight", provide a unique probe of the thermodynamics of these processes, with the capability to directly measure the density, pressure, and temperature of ionized gas. Moreover, these observations have uniquely high sensitivity into the outskirts of the halos of galaxies and clusters, including systems at high redshift. In the next decade, the combination of large spectroscopic and photometric optical galaxy surveys and wide-field, low-noise CMB surveys will transform our understanding of galaxy formation via these probes.

We introduce AREPO-RT, a novel radiation hydrodynamic (RHD) solver for the unstructured moving-mesh code AREPO. Our method solves the moment-based radiative transfer equations using the M1 closure relation. We achieve second order convergence by using a slope limited linear spatial extrapolation and a first order time prediction step to obtain the values of the primitive variables on both sides of the cell interface. A Harten-Lax-Van Leer flux function, suitably modified for moving meshes, is then used to solve the Riemann problem at the interface. The implementation is fully conservative and compatible with the individual timestepping scheme of AREPO. It incorporates atomic Hydrogen (H) and Helium (He) thermochemistry, which is used to couple the ultra-violet (UV) radiation field to the gas. Additionally, infrared radiation is coupled to the gas under the assumption of local thermodynamic equilibrium between the gas and the dust. We successfully apply our code to a large number of test problems, including applications such as the expansion of ${\rm H_{II}}$ regions, radiation pressure driven outflows and the levitation of optically thick layer of gas by trapped IR radiation. The new implementation is suitable for studying various important astrophysical phenomena, such as the effect of radiative feedback in driving galactic scale outflows, radiation driven dusty winds in high redshift quasars, or simulating the reionisation history of the Universe in a self consistent manner.

The "sloshing" of the cold gas in the cores of relaxed clusters of galaxies is a widespread phenomenon, evidenced by the presence of spiral-shaped "cold fronts" in X-ray observations of these systems. In simulations, these flows of cold gas readily form by interactions of the cluster core with small subclusters, due to a separation of the cold gas from the dark matter (DM), due to their markedly different collisionalities. In this work, we use numerical simulations to investigate the effects of increasing the DM collisionality on sloshing cold fronts in a cool-core cluster. For clusters in isolation, the formation of a flat DM core via self-interactions results in modest adiabatic expansion and cooling of the core gas. In merger simulations, cold fronts form in the same manner as in previous simulations, but the flattened potential in the core region enables the gas to expand to larger radii in the initial stages. Upon infall, the subcluster's DM mass decreases via collisions, reducing its influence on the core. Thus, the sloshing gas moves slower, inhibiting the growth of fluid instabilities relative to simulations where the DM cross section is zero. This also inhibits turbulent mixing and the increase in entropy that would otherwise result. For values of the cross section $\sigma/m > 1$, subclusters do not survive as self-gravitating structures for more than two core passages. Additionally, separations between the peaks in the X-ray emissivity and thermal Sunyaev-Zeldovich effect signals during sloshing may place constraints on DM self-interactions.

In this work, we discuss the relation between the strength of the self-interaction of dark matter particles and the predicted properties of the inner density distributions of dark matter haloes. We present the results of $N$-body simulations for 28 haloes performed with the same initial conditions for cold dark matter and for self-interacting dark matter, for a range of cross-sections. We provide a simple phenomenological description of these results and compare with the semi-analytical model typically used in the literature. Using these results, we then predict how the inner dark matter surface density and core radius should depend on the self-interaction cross-section for observed haloes.

We use high-resolution hydrodynamical simulation to test the difference of halo properties in cold dark matter (CDM) and a self-interacting dark matter (SIDM) scenario with a constant cross-section of $\sigma^\text{T}/m_{\chi}=1\;\text{cm}^{2}\text{g}^{-1}$. We find that the interplay between dark matter self-interaction and baryonic physics induces a complex evolution of the halo properties, which depends on the halo mass and morphological type, as well as on the halo mass accretion history. While high mass haloes, selected as analogues of early-type galaxies, show cored profiles in the SIDM run, systems of intermediate mass and with a significant disk component can develop a profile that is similar or cuspier than in CDM. The final properties of SIDM haloes - measured at z=0.2 - correlate with the halo concentration and formation time, suggesting that the differences between different systems are due to the fact that we are observing the impact self-interaction. We also search for signatures of self-interacting dark matter in the lensing signal of the main haloes and hints of potential differences in the distribution of Einstein radii, which suggests that future wide-field survey might be able to distinguish between CDM and SIDM models on this basis. Finally, we find that the subhalo abundances are not altered in the adopted SIDM model with respect to CDM.

We use the IllustrisTNG (TNG) cosmological simulations to provide theoretical expectations for the dark matter mass fractions (DMFs) and circular velocity profiles of galaxies. TNG predicts flat circular velocity curves for $z = 0$ Milky Way (MW)-like galaxies beyond a few kpc from the galaxy centre, in better agreement with observational constraints than its predecessor, Illustris. TNG also predicts an enhancement of the dark matter mass within the 3D stellar half-mass radius ($r_\mathrm{half}$; $M_\mathrm{200c} = 10^{10}-10^{13}\mathrm{M}_{\odot}$, $z \le2$) compared to its dark matter only and Illustris counterparts. This enhancement leads TNG present-day galaxies to be dominated by dark matter within their inner regions, with $f_\mathrm{DM}(<r_\mathrm{half})\gtrsim0.5$ at all masses and with a minimum for MW-mass galaxies. The 1$\sigma$ scatter is $\lesssim$ 10~per~cent at all apertures, which is smaller than that inferred by some observational datasets, e.g. 40 per cent from the SLUGGS survey. TNG agrees with the majority of the observationally inferred values for elliptical galaxies once a consistent IMF is adopted (Chabrier) and the DMFs are measured within the same apertures. The DMFs measured within $r_\mathrm{half}$ increase towards lower redshifts: this evolution is dominated by the increase in galaxy size with time. At $z\sim2$, the DMF in disc-like TNG galaxies decreases with increasing galaxy mass, with $f_\mathrm{DM}(<r_\mathrm{half}) \sim 0.10-0.65$ for $10^{10} \lesssim M_{\rm stars}/\mathrm{M}_{\odot} \lesssim 10^{12}$, and are two times higher than if TNG galaxies resided in Navarro-Frenk-White dark matter haloes unaffected by baryonic physics. It remains to be properly assessed whether recent observational estimates of the DMFs at $z\sim2$ rule out the contraction of the dark matter haloes predicted by the TNG model.

We contrast predictions for the high-redshift galaxy population and reionization history between cold dark matter (CDM) and an alternative self-interacting dark matter model based on the recently developed ETHOS framework that alleviates the small-scale CDM challenges within the Local Group. We perform the highest resolution hydrodynamical cosmological simulations (a 36~Mpc$^3$ volume with gas cell mass of $\sim10^5\mathrm{M}_{\odot}$ and minimum gas softening of $\sim180$~pc) within ETHOS to date -- plus a CDM counterpart -- to quantify the abundance of galaxies at high redshift and their impact on reionization. We find that ETHOS predicts galaxies with higher ultraviolet (UV) luminosities than their CDM counterparts and a faster build-up of the faint end of the UV luminosity function. These effects, however, make the optical depth to reionization less sensitive to the power spectrum cut-off: the ETHOS model differs from the CDM $\tau$ value by only 10 per cent and is consistent with Planck limits if the effective escape fraction of UV photons is 0.1-0.5. We conclude that current observations of high-redshift luminosity functions cannot differentiate between ETHOS and CDM models, but deep JWST surveys of strongly-lensed, inherently faint galaxies have the potential to test non-CDM models that offer attractive solutions to CDM's Local Group problems.

We study the connection between active galactic nuclei (AGN) and their host galaxies through cosmic time in the large-scale cosmological IllustrisTNG simulations. We first compare BH properties, i.e. the hard X-ray BH luminosity function, AGN galaxy occupation fraction, and distribution of Eddington ratios, to available observational constraints. The simulations produce a population of BHs in good agreement with observations, but we note an excess of faint AGN in hard X-ray (L_x ~ 10^{43-44} erg/s), and a lower number of bright AGN (L_x>10^{44} erg/s), a conclusion that varies quantitatively but not qualitatively with BH luminosity estimation method. The lower Eddington ratios of the 10^{9} Msun BHs compared to observations suggest that AGN feedback may be too efficient in this regime. We study galaxy star formation activity and structural properties, and design sample-dependent criteria to identify different galaxy types (star-forming/quiescent, extended/compact) that we apply both to the simulations and observations from the candels fields. We analyze how the simulated and observed galaxies populate the specific star formation rate - stellar mass surface density diagram. A large fraction of the z=0 M_{star}>10^{11} Msun quiescent galaxies first experienced a compaction phase (i.e. reduction of galaxy size) while still forming stars, and then a quenching event. We measure the dependence of AGN fraction on galaxies' locations in this diagram. After correcting the simulations with a redshift and AGN luminosity-dependent model for AGN obscuration, we find good qualitative and quantitative agreement with observations. The AGN fraction is the highest among compact star-forming galaxies (16-20% at z~1.5-2), and the lowest among compact quiescent galaxies (6-10% at z~1.5-2).

We have generated synthetic images of $\sim$27,000 galaxies from the IllustrisTNG and the original Illustris hydrodynamic cosmological simulations, designed to match Pan-STARRS observations of $\log_{10}(M_{\ast}/{\rm M}_{\odot}) \approx 9.8$-$11.3$ galaxies at $z \approx 0.05$. Most of our synthetic images were created with the SKIRT radiative transfer code, including the effects of dust attenuation and scattering, and performing the radiative transfer directly on the Voronoi mesh used by the simulations themselves. We have analysed both our synthetic and real Pan-STARRS images with the newly developed $\tt{statmorph}$ code, which calculates non-parametric morphological diagnostics -- including the Gini-$M_{20}$ and concentration-asymmetry-smoothness (CAS) statistics -- and performs two-dimensional S\'ersic fits. Overall, we find that the optical morphologies of IllustrisTNG galaxies are in good agreement with observations, and represent a substantial improvement compared to the original Illustris simulation. In particular, the locus of the Gini-$M_{20}$ diagram is consistent with that inferred from observations, while the median trends with stellar mass of all the morphological, size and shape parameters considered in this work lie within the $\sim$1$\sigma$ scatter of the observational trends. However, the IllustrisTNG model has some difficulty with more stringent tests, such as producing a strong morphology-colour relation. This results in a somewhat higher fraction of red discs and blue spheroids compared to observations. Similarly, the morphology-size relation is problematic: while observations show that discs tend to be larger than spheroids at a fixed stellar mass, such a trend is not present in IllustrisTNG.

We study the chaotic-like behavior of cosmological simulations by quantifying how minute perturbations grow over time and manifest as macroscopic differences in galaxy properties. When we run pairs of 'shadow' simulations that are identical except for random minute initial displacements to particle positions (e.g. of order 1e-7pc), the results diverge from each other at the individual galaxy level (while the statistical properties of the ensemble of galaxies are unchanged). After cosmological times, the global properties of pairs of 'shadow' galaxies that are matched between the simulations differ from each other generally at a level of ~2-25%, depending on the considered physical quantity. We perform these experiments using cosmological volumes of (25-50Mpc/h)^3 evolved either purely with dark matter, or with baryons and star-formation but no feedback, or using the full feedback model of the IllustrisTNG project. The runs cover four resolution levels spanning a factor of 512 in mass. We find that without feedback the differences between shadow galaxies generally become smaller as the resolution increases, but with the IllustrisTNG model the results are mostly converging towards a 'floor'. This hints at the role of feedback in setting the chaotic properties of galaxy formation. Importantly, we compare the macroscopic differences between shadow galaxies to the overall scatter in various galaxy scaling relations, and conclude that for the star formation-mass and the Tully-Fisher relations the butterfly effect in our simulations contributes significantly to the overall scatter. We find that our results are robust to whether random numbers are used in the sub-grid models or not. We discuss the implications for galaxy formation theory in general and for cosmological simulations in particular.

A cutoff in the linear matter power spectrum at dwarf galaxy scales has been shown to affect the abundance, formation mechanism and age of dwarf haloes and their galaxies at high and low redshift. We use hydrodynamical simulations of galaxy formation within the ETHOS framework in a benchmark model that has such a cutoff, and that has been shown to be an alternative to the cold dark matter (CDM) model that alleviates its dwarf-scale challenges. We show how galaxies in this model form differently to CDM on a halo-by-halo basis, at redshifts $z\ge6$. We show that ETHOS haloes at the half-mode mass scale form with 50~per~cent less mass than their CDM counterparts due to their later formation times, yet they retain more of their gas reservoir due to the different behaviour of gas and dark matter during the monolithic collapse of the first haloes in models with a galactic-scale cutoff. As a result, galaxies in ETHOS haloes near the cutoff scale grow rapidly between $z=10-6$ and by $z=6$ end up having very similar stellar masses, higher gas fractions and higher star formation rates relative to their CDM counterparts. We highlight these differences by making predictions for how the number of galaxies with old stellar populations is suppressed in ETHOS for both $z=6$ galaxies and for gas-poor Local Group fossil galaxies. Interestingly, we find an age gradient in ETHOS between galaxies that form in high and low density environments.

Feedback from accreting supermassive black holes, active galactic nuclei (AGN), is now a cornerstone of galaxy formation models. In this work, we present radiation-hydrodynamic simulations of radiative AGN feedback using the novel Arepo-RT. A central black hole emits radiation at a constant luminosity and drives an outflow via radiation pressure on dust grains. Utilising an isolated NFW halo we validate our setup in the single and multi-scattering regimes, with the simulated shock front propagation in excellent agreement with the expected analytic result. For a spherically symmetric NFW halo, an examination of the simulated outflow properties generated by radiative feedback demonstrates that they are lower than typically observed at a fixed AGN luminosity, regardless of the collimation of the radiation. We then explore the impact of a central disc galaxy and the assumed dust model on the outflow properties. The contraction of the halo during the galaxy's formation and modelling the production of dust grains results in a factor $100$ increase in the halo's optical depth. Radiation is then able to couple momentum more efficiently to the gas, driving a stronger shock and producing a mass-loaded $\sim10^{3}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$ outflow with a velocity of $\sim2000\,\mathrm{km}\,\mathrm{s}^{-1}$, in agreement with observations. However, the inclusion of dust destruction mechanisms, like thermal sputtering, leads to the rapid destruction of dust grains within the outflow, reducing its properties below typically observed values. We conclude that radiative AGN feedback can drive outflows, but a thorough numerical and physical treatment is required to assess its true impact.

We investigate the spatial distribution of Lyman alpha emitting galaxies (LAEs) at $z \approx 2.67$, selected from the NOAO Deep Wide-Field Survey (NDWFS), using two-point statistics and topological diagnostics adopted from network science. We measure the clustering length, $r_0 \approx 4 h^{-1}$ Mpc, and the bias, $b_{LAE} = 2.2^{+0.2}_{-0.1}$. Fitting the clustering with halo occupation distribution (HOD) models results in two disparate possibilities: (1) where the fraction of central galaxies is $<$1% in halos of mass $>10^{12}$$M_\odot$; and (2) where the fraction is $\approx$20%. We refer to these two scenarios as the `Dusty Core Scenario' for Model#1 since most of central galaxies in massive halos are dead in Ly$\alpha$ emission, and the `Pristine Core Scenario' for Model#2 since the central galaxies are bright in Ly$\alpha$ emission. Traditional two-point statistics cannot distinguish between these disparate models given the current data sets. To overcome this degeneracy, we generate mock catalogs for each HOD model using a high resolution $N$-body simulation and adopt a network statistics approach, which provides excellent topological diagnostics for galaxy point distributions. We find three topological anomalies from the spatial distribution of observed LAEs, which are not reproduced by the HOD mocks. We find that Model#2 matches better all network statistics than Model#1, suggesting that the central galaxies in $> 10^{12} h^{-1} M_\odot$ halos at $z \approx 2.67$ need to be less dusty to be bright as LAEs, potentially implying some replenishing channels of pristine gas such as the cold mode accretion.

We use IllustrisTNG, a suite of gravity and MHD simulations, to study the demographics and properties of jellyfish galaxies in the full cosmological context. By jellyfish galaxies, we mean satellites orbiting in massive groups and clusters that exhibit highly asymmetric distributions of gas and gas tails. We use the TNG100 run and select galaxies at redshifts $z\le0.6$ with stellar mass exceeding $10^{9.5}{\rm M_\odot}$ and with host halo masses of $10^{13}-10^{14.6}\,{\rm M_\odot}$. Among more than about 6000 (2600) galaxies with stars (and some gas), we identify 800 jellyfish galaxies by visually inspecting their gas and stellar mass maps in random projections. About $31\%$ of cluster satellites are found with signatures of ram-pressure stripping and gaseous tails stemming from the main luminous bodies. This is a lower limit, since the random orientation entails a loss of about $30\%$ of galaxies that in an optimal projection would otherwise be identified as jellyfish. The connection with ram-pressure stripping is further confirmed by a series of findings: jellyfish galaxies are more frequent at intermediate and large cluster-centric distances ($r/R_{\rm 200c}\gtrsim 0.25$); they move through the ICM with larger bulk velocities and Mach numbers than the general cluster population, typically orbiting supersonically and experiencing larger ram pressures. Furthermore, the gaseous tails usually extend in opposite directions to the galaxy trajectory, with no relation between tail orientation and the host's center. The frequency of jellyfish galaxies shows a very weak dependence on redshift $(0\le z\le0.6)$ but larger fractions of disturbed gaseous morphologies occur in more massive hosts and at smaller satellite masses. Finally, jellyfish galaxies are late infallers ($< 2.5-3$ Gyrs ago, at $z=0$) and the emergence of gaseous tails correlates well with the presence of bow shocks in the ICM.

The thermodynamic structure of hot gas in galaxy clusters is sensitive to astrophysical processes and typically difficult to model with galaxy formation simulations. We explore the fraction of cool-core (CC) clusters in a large sample of $370$ clusters from IllustrisTNG, examining six common CC definitions. IllustrisTNG produces continuous CC criteria distributions, the extremes of which are classified as CC and non-cool-core (NCC), and the criteria are increasingly correlated for more massive clusters. At $z=0$, the CC fractions for $2$ criteria are in reasonable agreement with the observed fractions but the other $4$ CC fractions are lower than observed. This result is partly driven by systematic differences between the simulated and observed gas fraction profiles. The simulated CC fractions with redshift show tentative agreement with the observed fractions, but linear fits demonstrate that the simulated evolution is steeper than observed. The conversion of CCs to NCCs appears to begin later and act more rapidly in the simulations. Examining the fraction of CCs and NCCs defined as relaxed we find no evidence that CCs are more relaxed, suggesting that mergers are not solely responsible for disrupting CCs. A comparison of the median thermodynamic profiles defined by different CC criteria shows that the extent to which they evolve in the cluster core is dependent on the CC criteria. We conclude that the thermodynamic structure of galaxy clusters in IllustrisTNG shares many similarities with observations, but achieving better agreement most likely requires modifications of the underlying galaxy formation model.

Galactic winds exhibit a multiphase structure that consists of hot-diffuse and cold-dense phases. Here we present high-resolution idealised simulations of the interaction of a hot supersonic wind with a cold cloud with the moving-mesh code arepo in setups with and without radiative cooling. We demonstrate that cooling causes clouds with sizes larger than the cooling length to fragment in two- and three-dimensional simulations (2D and 3D). We confirm earlier 2D simulations by McCourt et al. 2018 and highlight differences of the shattering processes of 3D clouds that are exposed to a hot wind. The fragmentation process is quantified with a friends-of-friends analysis of shattered cloudlets and density power spectra. Those show that radiative cooling causes the power spectral index to gradually increase when the initial cloud radius is larger than the cooling length and with increasing time until the cloud is fully dissolved in the hot wind. A resolution of around 1 pc is required to reveal the effect of cooling-induced fragmentation of a 100 pc outflowing cloud. Thus, state-of-the-art cosmological zoom simulations of the circumgalactic medium (CGM) fall short by orders of magnitudes from resolving this fragmentation process. This physics is, however, necessary to reliably model observed column densities and covering fractions of Lyman-$\alpha$ haloes, high-velocity clouds, and broad-line regions of active galactic nuclei.

Large-scale cosmological simulations of galaxy formation currently do not resolve the densities at which molecular hydrogen forms, implying that the atomic-to-molecular transition must be modeled either on the fly or in postprocessing. We present an improved postprocessing framework to estimate the abundance of atomic and molecular hydrogen and apply it to the IllustrisTNG simulations. We compare five different models for the atomic-to-molecular transition, including empirical, simulation-based, and theoretical prescriptions. Most of these models rely on the surface density of neutral hydrogen and the ultraviolet (UV) flux in the Lyman-Werner band as input parameters. Computing these quantities on the kiloparsec scales resolved by the simulations emerges as the main challenge. We show that the commonly used Jeans length approximation to the column density of a system can be biased and exhibits large cell-to-cell scatter. Instead, we propose to compute all surface quantities in face-on projections and perform the modeling in two dimensions. In general, the two methods agree on average, but their predictions diverge for individual galaxies and for models based on the observed midplane pressure of galaxies. We model the UV radiation from young stars by assuming a constant escape fraction and optically thin propagation throughout the galaxy. With these improvements, we find that the five models for the atomic-to-molecular transition roughly agree on average but that the details of the modeling matter for individual galaxies and the spatial distribution of molecular hydrogen. We emphasize that the estimated molecular fractions are approximate due to the significant systematic uncertainties.

We study the effect of baryonic processes on the shapes of dark matter (DM) haloes from Illustris, a suite of hydrodynamical (Illustris) and DM-only (Illustris-Dark) cosmological simulations performed with the moving-mesh code {\sc arepo}. DM halo shapes are determined using an iterative method based on the inertia tensor for a wide range of $z=0$ masses ($M_{200} = 1 \times 10^{11} - 3 \times 10^{14} M_\odot$). Convergence tests shows that the local DM shape profiles are converged only for $r > 9\epsilon$, $\epsilon$ being the Plummer-equivalent softening length, larger than expected. Haloes from non-radiative simulations (i.e. neglecting radiative processes, star formation, and feedback) exhibit no alteration in shapes from their DM-only counterparts: thus moving-mesh hydrodynamics alone is insufficient to cause differences in DM shapes. With the full galaxy-physics implementation, condensation of baryons results in significantly rounder and more oblate haloes, with the median minor-to-major axis ratio $\left< s\equiv c/a \right> \approx 0.7$, almost invariant throughout the halo and across halo masses. This somewhat improves the agreement between simulation predictions and observational estimates of the Milky Way halo shape. Consistently, the velocity anisotropy of DM is also reduced in Illustris, across halo masses and radii. Within the inner halo ($r=0.15 R_{200}$), both $s$ and $q$ (intermediate-to-major axis ratio) exhibit non-monotonicity with galaxy mass, peaking at $m_* \approx 10^{10.5-11} M_\odot$, which we find is due to the strong dependence of inner halo shape with galaxy formation efficiency. Baryons in Illustris affect the correlation of halo shape with halo properties, leading to a positive correlation of sphericity of MW-mass haloes with halo formation time and concentration, the latter being mildly more pronounced than in Illustris-Dark.

We introduce the "Supersonic Project," aimed at investigating the effects of the supersonic relative velocity between dark matter (DM) and baryons at high redshift using a combination of analytical calculations and cosmological simulations. In this paper, we study the effect of this stream velocity on the angular momentum of the first structures in the early Universe using simulations. We focus on DM haloes and their gas component as well as the recently predicted supersonically-induced gas objects (SIGOs) that arise as a result of the stream velocity phase shift. We find that the spin parameter of the gas component in these first haloes is increased with the stream velocity. Moreover, we find that when the stream velocity is taken into account, the angular momentum vectors of the DM component and the gas component are typically misaligned and this misalignment angle has a nearly isotropic distribution. The spin parameter value of the gas component is higher than in the no stream velocity case, which even in the absence of cooling, may result in more prolate objects. We also generalize the spin parameter to the SIGOs and find that they typically have a larger spin parameter with respect to their dark matter counterparts and that there is no correlation of the spin parameter and the prolateness of such structures. We speculate that SIGOs may be observed as very low luminosity objects in the early Universe and may serve as potential progenitors of Little Blue Dot-like systems.

Strong gravitational lensing has been identified as a promising astrophysical probe to study the particle nature of dark matter. In this paper we present a detailed study of the power spectrum of the projected mass density (convergence) field of substructure in a Milky Way-sized halo. This power spectrum has been suggested as a key observable that can be extracted from strongly lensed images and yield important clues about the matter distribution within the lens galaxy. We use two different $N$-body simulations from the ETHOS framework: one with cold dark matter and another with self-interacting dark matter and a cutoff in the initial power spectrum. Despite earlier works that identified $ k \gtrsim 100$ kpc$^{-1}$ as the most promising scales to learn about the particle nature of dark matter we find that even at lower wavenumbers - which are actually within reach of observations in the near future - we can gain important information about dark matter. Comparing the amplitude and slope of the power spectrum on scales $0.1 \lesssim k/$kpc$^{-1} \lesssim 10$ from lenses at different redshifts can help us distinguish between cold dark matter and other exotic dark matter scenarios that alter the abundance and central densities of subhalos. Furthermore, by considering the contribution of different mass bins to the power spectrum we find that subhalos in the mass range $10^7 - 10^8$ M$_{\odot}$ are on average the largest contributors to the power spectrum signal on scales $2 \lesssim k/$kpc$^{-1} \lesssim 15$, despite the numerous subhalos with masses $> 10^8$ M$_{\odot}$ in a typical lens galaxy. Finally, by comparing the power spectra obtained from the subhalo catalogs to those from the particle data in the simulation snapshots we find that the seemingly-too-simple halo model is in fact a fairly good approximation to the much more complex array of substructure in the lens.

We introduce the IllustrisTNG project, a new suite of cosmological magnetohydrodynamical simulations performed with the moving-mesh code AREPO employing an updated Illustris galaxy formation model. Here we focus on the general properties of magnetic fields and the diffuse radio emission in galaxy clusters. Magnetic fields are prevalent in galaxies, and their build-up is closely linked to structure formation. We find that structure formation amplifies the initial seed fields ($10^{-14}$ comoving Gauss) to the values observed in low-redshift galaxies ($1-10\,\mu{\rm G}$). The magnetic field topology is closely connected to galaxy morphology such that irregular fields are hosted by early-type galaxies, while large-scale, ordered fields are present in disc galaxies. Using two simple models for the energy distribution of relativistic electrons we predict the diffuse radio emission of $280$ clusters with a baryonic mass resolution of $1.1\times 10^{7}\,{\rm M_{\odot}}$, and generate mock observations for VLA, LOFAR, ASKAP and SKA. Our simulated clusters show extended radio emission, whose detectability correlates with their virial mass. We reproduce the observed scaling relations between total radio power and X-ray emission, $M_{500}$, and the Sunyaev-Zel'dovich $Y_{\rm 500}$ parameter. The radio emission surface brightness profiles of our most massive clusters are in reasonable agreement with VLA measurements of Coma and Perseus. Finally, we discuss the fraction of detected extended radio haloes as a function of virial mass and source count functions for different instruments. Overall our results agree encouragingly well with observations, but a refined analysis requires a more sophisticated treatment of relativistic particles in large-scale galaxy formation simulations.

We employ isolated N-body simulations to study the response of self-interacting dark matter (SIDM) halos in the presence of the baryonic potentials. Dark matter self-interactions lead to kinematic thermalization in the inner halo, resulting in a tight correlation between the dark matter and baryon distributions. A deep baryonic potential shortens the phase of SIDM core expansion and triggers core contraction. This effect can be further enhanced by a large self-scattering cross section. We find the final SIDM density profile is sensitive to the baryonic concentration and the strength of dark matter self-interactions. Assuming a spherical initial halo, we also study evolution of the SIDM halo shape together with the density profile. The halo shape at later epochs deviates from spherical symmetry due to the influence of the non-spherical disc potential, and its significance depends on the baryonic contribution to the total gravitational potential, relative to the dark matter one. In addition, we construct a multi-component model for the Milky Way, including an SIDM halo, a stellar disc and a bulge, and show it is consistent with observations from stellar kinematics and streams.

The galaxy Malin 1 contains the largest stellar disk known but the formation mechanism of this structure has been elusive. In this paper, we report a Malin 1 analogue in the 100 Mpc IllustrisTNG simulation and describe its formation history. At redshift zero, this massive galaxy, having a maximum circular velocity $V_{\rm max}$ of 430 ${\rm km\ s^{-1}}$, contains a 100 kpc gas/stellar disk with morphology similar to Malin 1. The simulated galaxy reproduces well many observed features of Malin 1's vast disk, including its stellar ages, metallicities, and gas rotation curve. We trace the extended disk back in time and find that a large fraction of the cold gas at redshift zero originated from the cooling of hot halo gas, triggered by the merger of a pair of intruding galaxies. Our finding provides a novel way to form large galaxy disks as extreme as Malin 1 within the current galaxy formation framework.

Understanding how baryonic processes shape the intracluster medium (ICM) is of critical importance to the next generation of galaxy cluster surveys. However, most models of structure formation neglect potentially important physical processes, like anisotropic thermal conduction (ATC). In this letter, we explore the impact of ATC on the prevalence of cool-cores (CCs) using 12 pairs of magnetohydrodynamical galaxy cluster simulations, simulated using the IllustrisTNG model with and without ATC. Although the impact of ATC varies from cluster to cluster and with CC criterion, its inclusion produces a systematic shift to larger CC fractions at z = 0 for all CC criteria considered. Additionally, the inclusion of ATC yields a flatter CC fraction redshift evolution, easing the tension with the observed evolution. With ATC included, the energy required for the central black hole to achieve self-regulation is reduced and the gas fraction in the cluster core increases, resulting in larger CC fractions. ATC makes the ICM unstable to perturbations and the increased efficiency of AGN feedback suggests that its inclusion results in a greater level of mixing in the ICM. Therefore, ATC is potentially an important physical process in reproducing the thermal structure of the ICM.

Interstellar dust is an important component of the galactic ecosystem, playing a key role in multiple galaxy formation processes. We present a novel numerical framework for the dynamics and size evolution of dust grains implemented in the moving-mesh hydrodynamics code AREPO suited for cosmological galaxy formation simulations. We employ a particle-based method for dust subject to dynamical forces including drag and gravity. The drag force is implemented using a second-order semi-implicit integrator and validated using several dust-hydrodynamical test problems. Each dust particle has a grain size distribution, describing the local abundance of grains of different sizes. The grain size distribution is discretised with a second-order piecewise linear method and evolves in time according to various dust physical processes, including accretion, sputtering, shattering, and coagulation. We present a novel scheme for stochastically forming dust during stellar evolution and new methods for sub-cycling of dust physics time-steps. Using this model, we simulate an isolated disc galaxy to study the impact of dust physical processes that shape the interstellar grain size distribution. We demonstrate, for example, how dust shattering shifts the grain size distribution to smaller sizes resulting in a significant rise of radiation extinction from optical to near-ultraviolet wavelengths. Our framework for simulating dust and gas mixtures can readily be extended to account for other dynamical processes relevant in galaxy formation, like magnetohydrodynamics, radiation pressure, and thermo-chemical processes.

[Abridged] We study the abundance and clustering properties of HI at redshifts $z\leqslant5$ using TNG100, a large state-of-the-art magneto-hydrodynamic simulation of a 75 Mpc/h box size. We show that most of the HI lies within dark matter halos and quantify the average HI mass hosted by halos of mass M at redshift z. We find that only halos with circular velocities larger than $\simeq$ 30 km/s contain HI. While the density profiles of HI exhibit a large halo-to-halo scatter, the mean profiles are universal across mass and redshift. The HI in low-mass halos is mostly located in the central galaxy, while in massive halos is concentrated in the satellites. We show that the HI and matter density probability distribution functions differ significantly. Our results point out that for small halos the HI bulk velocity goes in the same direction and has the same magnitude as the halo peculiar velocity, while in large halos differences show up. We find that halo HI velocity dispersion follows a power-law with halo mass. We find a complicated HI bias, with HI becoming non-linear already at $k=0.3$ h/Mpc at $z\gtrsim3$. Our simulation reproduces the DLAs bias value from observations. We find that the clustering of HI can be accurately reproduced by perturbative methods. We identify a new secondary bias, by showing that the clustering of halos depends not only on mass but also on HI content. We compute the amplitude of the HI shot-noise and find that it is small at all redshifts. We study the clustering of HI in redshift-space, and show that linear theory can explain the ratio between the monopoles in redshift- and real-space down to small scales at high redshift. We find that the amplitude of the Fingers-of-God effect is larger for HI than for matter. We point out that accurate 21 cm maps can be created from N-body or approximate simulations rather than full hydrodynamic simulations.

We explore the abundance, spatial distribution, and physical properties of the OVI, OVII, and OVIII ions of oxygen in circumgalactic and intergalactic media (the CGM, IGM, and WHIM). We use the TNG100 and TNG300 large volume cosmological magneto-hydrodynamical simulations. Modeling the ionization states of simulated oxygen, we find good agreement with observations of the low-redshift OVI column density distribution function (CDDF), and present its evolution for all three ions from z=0 to z=4. Producing mock quasar absorption line spectral surveys, we show that the IllustrisTNG simulations are fully consistent with constraints on the OVI content of the CGM from COS-Halos and other low redshift observations, producing columns as high as observed. We measure the total amount of mass and average column densities of each ion using hundreds of thousands of simulated galaxies spanning 10^11 < Mhalo/Msun < 10^15 corresponding to 10^9 < M*/Msun < 10^12 in stellar mass. Stacked radial profiles of OVI are computed in 3D number density and 2D projected column, decomposing into the 1-halo and 2-halo terms. Relating halo OVI to properties of the central galaxy, we find a correlation between the (g-r) color of a galaxy and the total amount of OVI in its CGM. In comparison to the COS-Halos finding, this leads to a dichotomy of columns around star-forming versus passive galaxies at fixed stellar (or halo) mass. We demonstrate that this correlation is a direct result of blackhole feedback associated with quenching and represents a causal consequence of galactic-scale baryonic feedback impacting the physical state of the circumgalactic medium.

The distribution of elements in galaxies provides a wealth of information about their production sites and their subsequent mixing into the interstellar medium. Here we investigate the distribution of elements within stars in the IllustrisTNG simulations. In particular, we analyze the abundance ratios of magnesium and europium in Milky Way-like galaxies from the TNG100 simulation (stellar masses ${\log} (M_\star / {\rm M}_\odot) \sim 9.7 - 11.2$). As abundances of magnesium and europium for individual stars in the Milky Way are observed across a variety of spatial locations and metallicities, comparison with the stellar abundances in our more than $850$ Milky Way-like galaxies provides stringent constraints on our chemical evolutionary methods. To this end we use the magnesium to iron ratio as a proxy for the effects of our SNII and SNIa metal return prescription, and a means to compare our simulated abundances to a wide variety of galactic observations. The europium to iron ratio tracks the rare ejecta from neutron star -- neutron star mergers, the assumed primary site of europium production in our models, which in turn is a sensitive probe of the effects of metal diffusion within the gas in our simulations. We find that europium abundances in Milky Way-like galaxies show no correlation with assembly history, present day galactic properties, and average galactic stellar population age. In general, we reproduce the europium to iron spread at low metallicities observed in the Milky Way, with the level of enhancement being sensitive to gas properties during redshifts $z \approx 2-4$. We show that while the overall normalization of [Eu/Fe] is susceptible to resolution and post-processing assumptions, the relatively large spread of [Eu/Fe] at low [Fe/H] when compared to that at high [Fe/H] is very robust.

In certain astrophysical systems the commonly employed ideal magnetohydrodynamics (MHD) approximation breaks down. Here, we introduce novel explicit and implicit numerical schemes of ohmic resistivity terms in the moving-mesh code AREPO. We include these non-ideal terms for two MHD techniques: the Powell 8-wave formalism and a constrained transport scheme, which evolves the cell-centred magnetic vector potential. We test our implementation against problems of increasing complexity, such as one- and two-dimensional diffusion problems, and the evolution of progressive and stationary Alfv\'en waves. On these test problems, our implementation recovers the analytic solutions to second-order accuracy. As first applications, we investigate the tearing instability in magnetized plasmas and the gravitational collapse of a rotating magnetized gas cloud. In both systems, resistivity plays a key role. In the former case, it allows for the development of the tearing instability through reconnection of the magnetic field lines. In the latter, the adopted (constant) value of ohmic resistivity has an impact on both the gas distribution around the emerging protostar and the mass loading of magnetically driven outflows. Our new non-ideal MHD implementation opens up the possibility to study magneto-hydrodynamical systems on a moving mesh beyond the ideal MHD approximation.

We use the IllustrisTNG simulations to investigate the evolution of the mass-metallicity relation (MZR) for star-forming cluster galaxies as a function of the formation history of their cluster host. The simulations predict an enhancement in the gas-phase metallicities of star-forming cluster galaxies (10^9< M_star<10^10 M_sun) at z<1.0 in comparisons to field galaxies. This is qualitatively consistent with observations. We find that the metallicity enhancement of cluster galaxies appears prior to their infall into the central cluster potential, indicating for the first time a systematic "chemical pre-processing" signature for {\it infalling} cluster galaxies. Namely, galaxies which will fall into a cluster by z=0 show a ~0.05 dex enhancement in the MZR compared to field galaxies at z<0.5. Based on the inflow rate of gas into cluster galaxies and its metallicity, we identify that the accretion of pre-enriched gas is the key driver of the chemical evolution of such galaxies, particularly in the stellar mass range (10^9< M_star<10^10 M_sun). We see signatures of an environmental dependence of the ambient/inflowing gas metallicity which extends well outside the nominal virial radius of clusters. Our results motivate future observations looking for pre-enrichment signatures in dense environments.

The IllustrisTNG project is a new suite of cosmological magneto-hydrodynamical simulations of galaxy formation performed with the Arepo code and updated models for feedback physics. Here we introduce the first two simulations of the series, TNG100 and TNG300, and quantify the stellar mass content of about 4000 massive galaxy groups and clusters ($10^{13} \leq M_{\rm 200c}/M_{\rm sun} \leq 10^{15}$) at recent times ($z \leq 1$). The richest clusters have half of their total stellar mass bound to satellite galaxies, with the other half being associated with the central galaxy and the diffuse intra-cluster light. The exact ICL fraction depends sensitively on the definition of a central galaxy's mass and varies in our most massive clusters between 20 to 40% of the total stellar mass. Haloes of $5\times 10^{14}M_{\rm sun}$ and above have more diffuse stellar mass outside 100 kpc than within 100 kpc, with power-law slopes of the radial mass density distribution as shallow as the dark matter's ( $-3.5 < \alpha_{\rm 3D} < -3$). Total halo mass is a very good predictor of stellar mass, and vice versa: at $z=0$, the 3D stellar mass measured within 30 kpc scales as $\propto (M_{\rm 500c})^{0.49}$ with a $\sim 0.12$ dex scatter. This is possibly too steep in comparison to the available observational constraints, even though the abundance of TNG less massive galaxies ($< 10^{11}M_{\rm sun}$ in stars) is in good agreement with the measured galaxy stellar mass functions at recent epochs. The 3D sizes of massive galaxies fall too on a tight ($\sim$0.16 dex scatter) power-law relation with halo mass, with $r^{\rm stars}_{\rm 0.5} \propto (M_{\rm 500c})^{0.53}$. Even more fundamentally, halo mass alone is a good predictor for the whole stellar mass profiles beyond the inner few kpc, and we show how on average these can be precisely recovered given a single mass measurement of the galaxy or its halo.

We introduce the first two simulations of the IllustrisTNG project, a next generation of cosmological magnetohydrodynamical simulations, focusing on the optical colors of galaxies. We explore TNG100, a rerun of the original Illustris box, and TNG300, which includes 2x2500^3 resolution elements in a volume twenty times larger. Here we present first results on the galaxy color bimodality at low redshift. Accounting for the attenuation of stellar light by dust, we compare the simulated (g-r) colors of 10^9 < M*/Msun < 10^12.5 galaxies to the observed distribution from the Sloan Digital Sky Survey (SDSS). We find a striking improvement with respect to the original Illustris simulation, as well as excellent quantitative agreement in comparison to the observations, with a sharp transition in median color from blue to red at a characteristic M* ~ 10^10.5 Msun. Investigating the build-up of the color-mass plane and the formation of the red sequence, we demonstrate that the primary driver of galaxy color transition in the TNG model is supermassive blackhole feedback in its low-accretion state. Across the entire population we measure a median color transition timescale dt_green of ~1.6 Gyr, a value which drops for increasingly massive galaxies. We find signatures of the physical process of quenching: at fixed stellar mass, redder galaxies have lower SFRs, gas fractions, and gas metallicities; their stellar populations are also older and their large-scale interstellar magnetic fields weaker than in bluer galaxies. Finally, we measure the amount of stellar mass growth on the red sequence. Galaxies with M* > 10^11 Msun which redden at z<1 accumulate on average ~25% of their final z=0 mass post-reddening; at the same time, ~18% of such massive galaxies acquire half or more of their final stellar mass while on the red sequence.

Modern cosmological simulations model the universe with increasing sophistication and at higher spatial and temporal resolutions. These enhancements permit detailed comparisons between the simulation outputs and real observational data. Recent projects such as Illustris are capable of producing simulated images that are comparable to those obtained from local surveys. This paper tests how well Illustris achieves this goal across a diverse population of galaxies using visual morphologies derived from Galaxy Zoo citizen scientists. Morphological classifications provided by volunteers for simulated galaxies are compared with similar data for a compatible sample of images drawn from the SDSS Legacy Survey. This paper investigates how simple morphological characterization by human volunteers asked to distinguish smooth from featured systems differs between simulated and real galaxy images. Differences are identified, which are likely due to the limited resolution of the simulation, but which could be revealing real differences in the dynamical evolution of populations of galaxies in the real and model universes. Specifically, for stellar masses $M_{\star}\lesssim10^{11}M_{\odot}$, a larger proportion of Illustris galaxies that exhibit disk-like morphology or visible substructure, relative to their SDSS counterparts. Toward higher masses, simulated and observed galaxies converge and exhibit similar morphology distributions. The stellar mass threshold indicated by this divergent behavior confirms recent works using parametric measures of morphology from Illustris simulated images. When $M_{\star}\gtrsim10^{11}M_{\odot}$, the Illustris dataset contains fewer galaxies that classifiers regard as unambiguously featured. These results suggest that comparison between the detailed properties of observed and simulated galaxies, even when limited to reasonably massive systems, may be misleading.

Galaxies are surrounded by sizeable gas reservoirs which host a significant amount of metals: the circum-galactic medium (CGM). The CGM acts as a mediator between the galaxy and the extra-galactic medium. However, our understanding of how galaxy mergers, a major evolutionary transformation, impact the CGM remains deficient. We present a theoretical study of the effect of galaxy mergers on the CGM. We use hydrodynamical cosmological zoom-in simulations of a major merger selected from the Illustris project such that the z=0 descendant has a halo mass and stellar mass comparable to the Milky Way. To study the CGM we then re-simulated this system at a 40 times better mass resolution, and included detailed post-processing ionization modelling. Our work demonstrates the effect the merger has on the characteristic size of the CGM, its metallicity, and the predicted covering fraction of various commonly observed gas-phase species, such as H I, C IV, and O VI. We show that merger-induced outflows can increase the CGM metallicity by 0.2-0.3 dex within 0.5 Gyr post-merger. These effects last up to 6 Gyr post-merger. While the merger increases the total metal covering fractions by factors of 2-3, the covering fractions of commonly observed UV ions decrease due to the hard ionizing radiation from the active galactic nucleus, which we model explicitly. Our study of the single simulated major merger presented in this work demonstrates the significant impact that a galaxy interaction can have on the size, metallicity, and observed column densities of the CGM.

Hydrodynamical simulations of galaxy formation have now reached sufficient volume to make precision predictions for clustering on cosmologically relevant scales. Here we use our new IllustrisTNG simulations to study the non-linear correlation functions and power spectra of baryons, dark matter, galaxies and haloes over an exceptionally large range of scales. We find that baryonic effects increase the clustering of dark matter on small scales and damp the total matter power spectrum on scales up to k ~ 10 h/Mpc by 20%. The non-linear two-point correlation function of the stellar mass is close to a power-law over a wide range of scales and approximately invariant in time from very high redshift to the present. The two-point correlation function of the simulated galaxies agrees well with SDSS at its mean redshift z ~ 0.1, both as a function of stellar mass and when split according to galaxy colour, apart from a mild excess in the clustering of red galaxies in the stellar mass range 10^9-10^10 Msun/h^2. Given this agreement, the TNG simulations can make valuable theoretical predictions for the clustering bias of different galaxy samples. We find that the clustering length of the galaxy auto-correlation function depends strongly on stellar mass and redshift. Its power-law slope gamma is nearly invariant with stellar mass, but declines from gamma ~ 1.8 at redshift z=0 to gamma ~ 1.6 at redshift z ~ 1, beyond which the slope steepens again. We detect significant scale-dependencies in the bias of different observational tracers of large-scale structure, extending well into the range of the baryonic acoustic oscillations and causing nominal (yet fortunately correctable) shifts of the acoustic peaks of around ~5%.

The fundamental metallicity relation (FMR) is a postulated correlation between galaxy stellar mass, star formation rate (SFR), and gas-phase metallicity. At its core, this relation posits that offsets from the mass-metallicity relation (MZR) at a fixed stellar mass are correlated with galactic SFR. In this Letter, we quantify the timescale with which galactic SFRs and metallicities evolve using hydrodynamical simulations. We find that Illustris and IllustrisTNG predict that galaxy offsets from the star formation main sequence and MZR evolve over similar timescales, are often anti-correlated in their evolution, evolve with the halo dynamical time, and produce a pronounced FMR. In fact, for a FMR to exist, the metallicity and SFR must evolve in an anti-correlated sense which requires that they evolve with similar time variability. In contrast to Illustris and IllustrisTNG, we speculate that the SFR and metallicity evolution tracks may become decoupled in galaxy formation models dominated by globally-bursty SFR histories, which could weaken the FMR residual correlation strength. This opens the possibility of discriminating between bursty and non-bursty feedback models based on the strength and persistence of the FMR -- especially at high redshift.

We analyze scaling relations and evolution histories of galaxy sizes in TNG100, part of the IllustrisTNG simulation suite. Observational qualitative trends of size with stellar mass, star-formation rate and redshift are reproduced, and a quantitative comparison of projected r-band sizes at 0~<z<~2 shows agreement to much better than 0.25dex. We follow populations of z=0 galaxies with a range of masses backwards in time along their main progenitor branches, distinguishing between main-sequence and quenched galaxies. Our main findings are as follows. (i) At M_{*,z=0}>~10^{9.5}Msun, the evolution of the median main progenitor differs, with quenched galaxies hardly growing in median size before quenching, whereas main-sequence galaxies grow their median size continuously, thus opening a gap from the progenitors of quenched galaxies. This is partly because the main-sequence high-redshift progenitors of quenched z=0 galaxies are drawn from the lower end of the size distribution of the overall population of main-sequence high-redshift galaxies. (ii) Quenched galaxies with M_{*,z=0}>~10^{9.5}Msun experience a steep size growth on the size-mass plane after their quenching time, but with the exception of galaxies with M_{*,z=0}>~10^{11}Msun, the size growth after quenching is small in absolute terms, such that most of the size (and mass) growth of quenched galaxies (and its variation among them) occurs while they are still on the main-sequence. After they become quenched, the size growth rate of quenched galaxies as a function of time, as opposed to versus mass, is similar to that of main-sequence galaxies. Hence, the size gap is retained down to z=0.

X-ray emitting gaseous coronae around massive galaxies are a basic prediction of galaxy formation models. Although the coronae around spiral galaxies offer a fundamental test of these models, observational constraints on their characteristics are still scarce. While the presence of extended hot coronae has been established around a handful of massive spiral galaxies, the short X-ray observations only allowed for measurements of the basic characteristics of the coronae. In this work, we utilize deep XMM-Newton observations of NGC 6753 to explore its extended X-ray corona in unprecedented detail. Specifically, we establish the isotropic morphology of the hot gas, suggesting that it resides in hydrostatic equilibrium. The temperature profile of the gas shows a decrease with increasing radius: it drops from $kT\approx0.7$ keV in the innermost parts to $kT\approx0.4$ keV at 50 kpc radius. The temperature map reveals the complex temperature structure of the gas. We study the metallicity distribution of the gas, which is uniform at $Z\approx0.1$ Solar. This value is about an order of magnitude lower than that obtained for elliptical galaxies with similar dark matter halo mass, hinting that the hot gas in spiral galaxies predominantly originates from external gas inflows rather than from internal sources. By extrapolating the density profile of the hot gas out to the virial radius, we estimate the total gas mass and derive the total baryon mass of NGC 6753. We conclude that the baryon mass fraction is $f_{\rm b} \approx 0.06$, implying that about half of the baryons are missing.

The coevolution of galaxies and their metal content serves as an important test for galaxy feedback models. We analyze the distribution and evolution of metals within the IllustrisTNG simulation suite with a focus on the gas-phase mass-metallicity relation (MZR). We find that the IllustrisTNG model broadly reproduces the slope and normalization evolution of the MZR across the redshift range $0<z<2$ and mass range $10^9 < M_*/\mathrm{M}_\odot < 10^{10.5}$. We make predictions for the high redshift ($2<z<10$) metal content of galaxies which is described by a gradual decline in the normalization of the metallicity with an average high redshift ($z>2$) evolution fit by $\mathrm{d\;log(Z)}/\mathrm{dz} \approx - 0.064$. Our simulations indicate that the metal retention efficiency of the interstellar medium (ISM) is low: a majority of gas-phase metals ($\sim$ 85 per cent at $z=0$) live outside of the ISM, either in an extended gas disk, the circumgalactic medium, or outside the halo. Nevertheless, the redshift evolution in the simulated MZR normalization is driven by the higher gas fractions of high redshift galaxies, not by changes to the metal retention efficiency. The scatter in the simulated MZR contains a clear correlation with the gas-mass or star formation rate of the system, in agreement with the observed fundamental metallicity relation. The scatter in the MZR is driven by a competition between periods of enrichment- and accretion-dominated metallicity evolution. We expect that while the normalization of the MZR declines with redshift, the strength of the correlation between metallicity and gas-mass at fixed stellar mass is not a strong function of redshift. Our results indicate that the "regulator" style models are best suited for simultaneously explaining the shape, redshift evolution, and existence of correlated scatter with gas fraction about the MZR.

The distribution of metals in the intra-cluster medium encodes important information about the enrichment history and formation of galaxy clusters. Here we explore the metal content of clusters in IllustrisTNG - a new suite of galaxy formation simulations building on the Illustris project. Our cluster sample contains 20 objects in TNG100 - a ~(100 Mpc)^3 volume simulation with 2x1820^3 resolution elements, and 370 objects in TNG300 - a ~(300 Mpc)^3 volume simulation with 2x2500^3 resolution elements. The z=0 metallicity profiles agree with observations, and the enrichment history is consistent with observational data going beyond z~1, showing nearly no metallicity evolution. The abundance profiles vary only minimally within the cluster samples, especially in the outskirts with a relative scatter of ~15%. The average metallicity profile flattens towards the center, where we find a logarithmic slope of -0.1 compared to -0.5 in the outskirts. Cool core clusters have more centrally peaked metallicity profiles (~0.8 solar) compared to non-cool core systems (~0.5 solar), similar to observational trends. Si/Fe and O/Fe radial profiles follow positive gradients. The outer abundance profiles do not evolve below z~2, whereas the inner profiles flatten towards z=0. More than ~80% of the metals in the intra-cluster medium have been accreted from the proto-cluster environment, which has been enriched to ~0.1 solar already at z~2. We conclude that the intra-cluster metal distribution is uniform among our cluster sample, nearly time-invariant in the outskirts for more than 10 Gyr, and forms through a universal enrichment history.

We introduce an updated physical model to simulate the formation and evolution of galaxies in cosmological, large-scale gravity+magnetohydrodynamical simulations with the moving mesh code AREPO. The overall framework builds upon the successes of the Illustris galaxy formation model, and includes prescriptions for star formation, stellar evolution, chemical enrichment, primordial and metal-line cooling of the gas, stellar feedback with galactic outflows, and black hole formation, growth and multi-mode feedback. In this paper we give a comprehensive description of the physical and numerical advances which form the core of the IllustrisTNG (The Next Generation) framework. We focus on the revised implementation of the galactic winds, of which we modify the directionality, velocity, thermal content, and energy scalings, and explore its effects on the galaxy population. As described in earlier works, the model also includes a new black hole driven kinetic feedback at low accretion rates, magnetohydrodynamics, and improvements to the numerical scheme. Using a suite of (25 Mpc $h^{-1}$)$^3$ cosmological boxes we assess the outcome of the new model at our fiducial resolution. The presence of a self-consistently amplified magnetic field is shown to have an important impact on the stellar content of $10^{12} M_{\rm sun}$ haloes and above. Finally, we demonstrate that the new galactic winds promise to solve key problems identified in Illustris in matching observational constraints and affecting the stellar content and sizes of the low mass end of the galaxy population.

We perform dark-matter-only simulations of 28 relaxed massive cluster-sized haloes for Cold Dark Matter (CDM) and Self-Interacting Dark Matter (SIDM) models, to study structural differences between the models at large radii, where the impact of baryonic physics is expected to be very limited. We find that the distributions for the radial profiles of the density, ellipsoidal axis ratios, and velocity anisotropies ($\beta$) of the haloes differ considerably between the models (at the $\sim1\sigma$ level), even at $\gtrsim10\%$ of the virial radius, if the self-scattering cross section is $\sigma/m_\chi=1$ cm$^2$ gr$^{-1}$. Direct comparison with observationally inferred density profiles disfavours SIDM for $\sigma/m_\chi=1$ cm$^2$ gr$^{-1}$, but in an intermediate radial range ($\sim3\%$ of the virial radius), where the impact of baryonic physics is uncertain. At this level of the cross section, we find a narrower $\beta$ distribution in SIDM, clearly skewed towards isotropic orbits, with no SIDM (90\% of CDM) haloes having $\beta>0.12$ at $7\%$ of the virial radius. We estimate that with an observational sample of $\sim30$ ($\sim10^{15}$ M$_\odot$) relaxed clusters, $\beta$ can potentially be used to put competitive constraints on SIDM, once observational uncertainties improve by a factor of a few. We study the suppression of the memory of halo assembly history in SIDM clusters. For $\sigma/m_\chi=1$ cm$^2$ gr$^{-1}$, we find that this happens only in the central halo regions ($\sim1/4$ of the scale radius of the halo), and only for haloes that assembled their mass within this region earlier than a formation redshift $z_f\sim2$. Otherwise, the memory of assembly remains and is reflected in ways similar to CDM, albeit with weaker trends.

We study the population of supermassive black holes (SMBHs) and their effects on massive central galaxies in the IllustrisTNG cosmological hydrodynamical simulations of galaxy formation. The employed model for SMBH growth and feedback assumes a two-mode scenario in which the feedback from active galactic nuclei occurs through a kinetic, comparatively efficient mode at low accretion rates relative to the Eddington limit, and in the form of a thermal, less efficient mode at high accretion rates. We show that the quenching of massive central galaxies happens coincidently with kinetic-mode feedback, consistent with the notion that active supermassive black cause the low specific star formation rates observed in massive galaxies. However, major galaxy mergers are not responsible for initiating most of the quenching events in our model. Up to black hole masses of about $10^{8.5}\,{\rm M}_\odot$, the dominant growth channel for SMBHs is in the thermal mode. Higher mass black holes stay mainly in the kinetic mode and gas accretion is self-regulated via their feedback, which causes their Eddington ratios to drop, with SMBH mergers becoming the main channel for residual mass growth. As a consequence, the quasar luminosity function is dominated by rapidly accreting, moderately massive black holes in the thermal mode. We show that the associated growth history of SMBHs produces a low-redshift quasar luminosity function and a redshift zero black hole mass-stellar bulge mass relation in good agreement with observations, whereas the simulation tends to over-predict the high-redshift quasar luminosity function.

We formulate an effective theory of structure formation (ETHOS) that enables cosmological structure formation to be computed in almost any microphysical model of dark matter physics. This framework maps the detailed microphysical theories of particle dark matter interactions into the physical effective parameters that shape the linear matter power spectrum and the self-interaction transfer cross section of non-relativistic dark matter. These are the input to structure formation simulations, which follow the evolution of the cosmological and galactic dark matter distributions. Models with similar effective parameters in ETHOS but with different dark particle physics would nevertheless result in similar dark matter distributions. We present a general method to map an ultraviolet complete or effective field theory of low energy dark matter physics into parameters that affect the linear matter power spectrum and carry out this mapping for several representative particle models. We further propose a simple but useful choice for characterizing the dark matter self-interaction transfer cross section that parametrizes self-scattering in structure formation simulations. Taken together, these effective parameters in ETHOS allow the classification of dark matter theories according to their structure formation properties rather than their intrinsic particle properties, paving the way for future simulations to span the space of viable dark matter physics relevant for structure formation.

We present a theoretical analysis of some unexplored aspects of relaxed Bose-Einstein condensate dark matter (BECDM) haloes. This type of ultralight bosonic scalar field dark matter is a viable alternative to the standard cold dark matter (CDM) paradigm, as it makes the same large-scale predictions as CDM and potentially overcomes CDM's small-scale problems via a galaxy-scale de Broglie wavelength. We simulate BECDM halo formation through mergers, evolved under the Schr\"odinger-Poisson equations. The formed haloes consist of a soliton core supported against gravitational collapse by the quantum pressure tensor and an asymptotic $r^{-3}$ NFW-like profile. We find a fundamental relation of the core=to-halo mass with the dimensionless invariant $\Xi \equiv \lvert E \rvert/M^3/(Gm/\hbar)^2$ or $M_{\rm c}/M \simeq 2.6 \Xi^{1/3}$, linking the soliton to global halo properties. For $r \geq 3.5 \,r_{\rm c}$ core radii, we find equipartition between potential, classical kinetic, and quantum gradient energies. The haloes also exhibit a conspicuous turbulent behavior driven by the continuous reconnection of vortex lines due to wave interference. We analyse the turbulence 1D velocity power spectrum and find a $k^{-1.1}$ power-law. This suggests the vorticity in BECDM haloes is homogeneous, similar to thermally-driven counterflow BEC systems from condensed matter physics, in contrast to a $k^{-5/3}$ Kolmogorov power-law seen in mechanically-driven quantum systems. The mode where the power spectrum peaks is approximately the soliton width, implying the soliton-sized granules carry most of the turbulent energy in BECDM haloes.

We present full volume cosmological simulations using the moving-mesh code AREPO to study the coevolution of dust and galaxies. We extend the dust model in AREPO to include thermal sputtering of grains and investigate the evolution of the dust mass function, the cosmic distribution of dust beyond the interstellar medium, and the dependence of dust-to-stellar mass ratio on galactic properties. The simulated dust mass function is well-described by a Schechter fit and lies closest to observations at $z = 0$. The radial scaling of projected dust surface density out to distances of $10 \, \text{Mpc}$ around galaxies with magnitudes $17 < i < 21$ is similar to that seen in Sloan Digital Sky Survey data, albeit with a lower normalisation. At $z = 0$, the predicted dust density of $\Omega_\text{dust} \approx 1.3 \times 10^{-6}$ lies in the range of $\Omega_\text{dust}$ values seen in low-redshift observations. We find that dust-to-stellar mass ratio anti-correlates with stellar mass for galaxies living along the star formation main sequence. Moreover, we estimate the $850 \, \mu\text{m}$ number density functions for simulated galaxies and analyse the relation between dust-to-stellar flux and mass ratios at $z = 0$. At high redshift, our model fails to produce enough dust-rich galaxies, and this tension is not alleviated by adopting a top-heavy initial mass function. We do not capture a decline in $\Omega_\text{dust}$ from $z = 2$ to $z = 0$, which suggests that dust production mechanisms more strongly dependent on star formation may help to produce the observed number of dusty galaxies near the peak of cosmic star formation.

We explore for the first time the effect of self-interacting dark matter (SIDM) on the dark matter (DM) and baryonic distribution in massive galaxies formed in hydrodynamical cosmological simulations, including explicit baryonic physics treatment. A novel implementation of Super-Massive Black Hole (SMBH) formation and evolution is used, as in Tremmel et al.(2015, 2016), allowing to explicitly follow SMBH dynamics at the center of galaxies. A high SIDM constant cross-section is chosen, $\sigma$=10 $\rm cm^2/gr$, to amplify differences from CDM models. Milky Way-like galaxies form a shallower DM density profile in SIDM than they do in CDM, with differences already at 20 kpc scales. This demonstrates that even for the most massive spirals the effect of SIDM dominates over the adiabatic contraction due to baryons. Strikingly, the dynamics of SMBHs differs in the SIDM and reference CDM case. SMBHs in massive spirals have sunk to the centre of their host galaxy in both the SIDM and CDM run, while in less massive galaxies about 80$\%$ of the SMBH population is off-centered in the SIDM case, as opposed to the CDM case in which $\sim$90$\%$ of SMBHs have reached their host's centre. SMBHs are found as far as $\sim$9 kpc away from the centre of their host SIDM galaxy. This difference is due to the increased dynamical friction timescale caused by the lower DM density in SIDM galaxies compared to CDM, resulting in 'core stalling'. This pilot work highlights the importance of simulating in a full hydrodynamical context different DM models combined to SMBH physics to study their influence on galaxy formation.

Early-type galaxies provide unique tests for the predictions of the cold dark matter cosmology and the baryonic physics assumptions entering models for galaxy formation. In this work, we use the Illustris simulation to study correlations of three main properties of early-type galaxies, namely, the stellar orbital anisotropies, the central dark matter fractions and the central radial density slopes, as well as their redshift evolution since $z=1.0$. We find that lower-mass galaxies or galaxies at higher redshift tend to be bluer in rest-frame colour, have higher central gas fractions, and feature more tangentially anisotropic orbits and steeper central density slopes than their higher-mass or lower-redshift counterparts, respectively. The projected central dark matter fraction within the effective radius shows a very mild mass dependence but positively correlates with galaxy effective radii due to the aperture effect. The central density slopes obtained by combining strong lensing measurements with single aperture kinematics are found to differ from the true density slopes. We identify systematic biases in this measurement to be due to two common modelling assumptions, isotropic stellar orbital distributions and power-law density profiles. We also compare the properties of early-type galaxies in Illustris to those from existing galaxy and strong lensing surveys, we find in general broad agreement but also some tension, which poses a potential challenge to the stellar formation and feedback models adopted by the simulation.

A power-law density model, i.e., $\rho(r) \propto r^{-\gamma'}$ has been commonly employed in strong gravitational lensing studies, including the so-called time-delay technique used to infer the Hubble constant $H_0$. However, since the radial scale at which strong lensing features are formed corresponds to the transition from the dominance of baryonic matter to dark matter, there is no known reason why galaxies should follow a power law in density. The assumption of a power law artificially breaks the mass-sheet degeneracy, a well-known invariance transformation in gravitational lensing which affects the product of Hubble constant and time delay and can therefore cause a bias in the determination of $H_0$ from the time-delay technique. In this paper, we use the Illustris hydrodynamical simulations to estimate the amplitude of this bias, and to understand how it is related to observational properties of galaxies. Investigating a large sample of Illustris galaxies that have velocity dispersion $\sigma_{SIE}$>160 km/s at redshifts below $z=1$, we find that the bias on $H_0$ introduced by the power-law assumption can reach 20%-50%, with a scatter of $10\%-30\%$ (rms). However, we find that by selecting galaxies with an inferred power-law model slope close to isothermal, it is possible to reduce the bias on $H_0$ to <5%, and the scatter to <10%. This could potentially be used to form less biased statistical samples for $H_0$ measurements in the upcoming large survey era.

Galactic rotation curves are a fundamental constraint for any cosmological model. We use controlled N-body simulations of galaxies to study the gravitational effect of baryons in a scenario with collisionless cold dark matter (CDM) versus one with a self-interacting dark matter (SIDM) component. In particular, we examine the inner profiles of the rotation curves in the velocity range Vmax = 30-250 km/s, whose diversity has been found to be greater than predicted by the Lambda-CDM scenario. We find that the scatter in the observed rotation curves exceeds that predicted by dark matter only mass-concentration relations in either the CDM nor SIDM models. Allowing for realistic baryonic content and spatial distributions, however, helps create a large variety of rotation curve shapes; which is in better agreement with observations in the case of self-interactions due to the characteristic cored profiles being more accommodating to the slowly rising rotation curves than CDM. We find individual fits to model two of the most remarkable outliers of similar Vmax, UGC 5721 and IC 2574, the former a cusp-like rotation curve and the latter a seemingly 8 kpc cored profile. This diversity in SIDM arises as permutations of overly concentrated halos with compact baryonic distributions versus underdense halos with extended baryonic disks. The SIDM solution is promising and its feasibility ultimately depends on the sampling of the halo mass-concentration relation and its interplay with the baryonic profiles, emphasising the need for a better understanding of the frequency of extreme outliers present in current observational samples.

We map the lensing-inferred substructure in the first three clusters observed by the Hubble Space Telescope Frontier Fields Initiative (HSTFF): Abell 2744 (z = 0.308), MACSJ0416, (z = 0.396) and MACSJ1149 (z = 0.543). Statistically resolving dark-matter subhaloes down to ~10^{9.5} solar masses, we compare the derived subhalo mass functions (SHMFs) to theoretical predictions from analytical models and with numerical simulations in a Lambda Cold Dark Matter (LCDM) cosmology. Mimicking our observational cluster member selection criteria in the HSTFF, we report excellent agreement in both amplitude and shape of the SHMF over four decades in subhalo mass (10^{9-13} solar masses). Projection effects do not appear to introduce significant errors in the determination of SHMFs from simulations. We do not find evidence for a substructure crisis, analogous to the missing satellite problem in the Local Group, on cluster scales, but rather excellent agreement of the count-matched HSTFF SHMF down to M_{sub halo}/M_{halo} ~ 10^{-5}. However, we do find discrepancies in the radial distribution of sub haloes inferred from HSTFF cluster lenses compared to determinations from simulated clusters. This suggests that although the selected simulated clusters match the HSTFF sample in mass, they do not adequately capture the dynamical properties and complex merging morphologies of these observed cluster lenses. Therefore, HSTFF clusters are likely observed in a transient evolutionary stage that is presently insufficiently sampled in cosmological simulations. The abundance and mass function of dark matter substructure in cluster lenses continues to offer an important test of the LCDM paradigm, and at present we find no tension between model predictions and observations.

We study intrinsic alignments (IA) of galaxy image shapes within the Illustris cosmic structure formation simulations. We investigate how IA correlations depend on observable galaxy properties such as stellar mass, apparent magnitude, redshift, and photometric type, and on the employed shape measurement method. The correlations considered include the matter density-intrinsic ellipticity (mI), galaxy density-intrinsic ellipticity (dI), gravitational shear-intrinsic ellipticity (GI), and intrinsic ellipticity-intrinsic ellipticity (II) correlations. We find stronger correlations for more massive and more luminous galaxies, as well as for earlier photometric types, in agreement with observations. Moreover, the correlations significantly depend on the choice of shape estimator, even if calibrated to serve as unbiased shear estimators. In particular, shape estimators that down-weight the outer parts of galaxy images produce much weaker IA signals on intermediate and large scales than methods employing flat radial weights. The expected contribution of intrinsic alignments to the observed ellipticity correlation in tomographic cosmic shear surveys may be below one percent or several percent of the full signal depending on the details of the shape measurement method. A comparison of our results to a tidal alignment model indicates that such a model is able to reproduce the IA correlations well on intermediate and large scales, provided the effect of varying galaxy density is correctly taken into account. We also find that the GI contributions to the observed ellipticity correlations could be inferred directly from measurements of galaxy density-intrinsic ellipticity correlations, except on small scales, where systematic differences between mI and dI correlations are large.

The inefficiency of star formation in massive elliptical galaxies is widely believed to be caused by the interactions of an active galactic nucleus (AGN) with the surrounding gas. Achieving a sufficiently rapid reddening of moderately massive galaxies without expelling too many baryons has however proven difficult for hydrodynamical simulations of galaxy formation, prompting us to explore a new model for the accretion and feedback effects of supermassive black holes. For high accretion rates relative to the Eddington limit, we assume that a fraction of the accreted rest mass energy heats the surrounding gas thermally, similar to the `quasar mode' in previous work. For low accretion rates, we invoke a new, pure kinetic feedback model which imparts momentum into the surrounding gas in a stochastic manner. These two modes of feedback are motivated both by theoretical conjectures for the existence of different types of accretion flows as well as recent observational evidence for the importance of kinetic AGN winds in quenching galaxies. We find that a large fraction of the injected kinetic energy in this mode thermalises via shocks in the surrounding gas, thereby providing a distributed heating channel. In cosmological simulations, the resulting model produces red, non star-forming massive elliptical galaxies, and achieves realistic gas fractions, black hole growth histories and thermodynamic profiles in large haloes.

We confirm that the object DDO216-A1 is a substantial globular cluster at the center of Local Group galaxy DDO216 (the Pegasus dwarf irregular), using Hubble Space Telescope ACS imaging. By fitting isochrones, we find the cluster metallicity to be [M/H] = -1.6 +/-0.2, for reddening E(B-V) = 0.16 +/-0.02; the best-fit age is 12.3 +/-0.8 Gyr. There are ~30 RR Lyrae variables in the cluster; the magnitude of the fundamental mode pulsators gives a distance modulus of 24.77 +/-0.08 - identical to the host galaxy. The ratio of overtone to fundamental mode variables and their mean periods make DDO216-A1 an Oosterhoff Type I cluster. We find an I-band central surface brightness 20.85 +/-0.17 F814W mag per square arcsecond, a half-light radius of 3.1 arcsec (13.4 pc), and an absolute magnitude M814 = -7.90 +/-0.16 (approximately 10^5 solar masses). King models fit to the cluster give the core radius and concentration index, r_c = 2.1" +/-0.9" and c = 1.24 +/-0.39. The cluster is an "extended" cluster somewhat typical of some dwarf galaxies and the outer halo of the Milky Way. The cluster is projected <30 pc south of the center of DDO216, unusually central compared to most dwarf galaxy globular clusters. Analytical models of dynamical friction and tidal destruction suggest that it probably formed at a larger distance, up to ~1 kpc, and migrated inward. DDO216 has an unexceptional cluster specific frequency, S_N = 10. DDO216 is the lowest-luminosity Local Group galaxy to host a 10^5 solar mass globular cluster, and the only transition-type (dSph/dIrr) in the Local Group with a globular.

Mergers and the spin of the dark matter halo are factors traditionally believed to determine the morphology of galaxies within a $\Lambda$CDM cosmology. We study this hypothesis by considering approximately 18,000 central galaxies at $z=0$ with stellar masses $M_{\ast} = 10^{9}-10^{12} \, {\rm M}_{\odot}$ selected from the Illustris cosmological hydrodynamic simulation. The fraction of accreted stars -- which measures the importance of massive, recent and dry mergers -- increases steeply with galaxy stellar mass, from less than 5 per cent in dwarfs to 80 per cent in the most massive objects, and the impact of mergers on galaxy morphology increases accordingly. For galaxies with $M_{\ast} \gtrsim 10^{11} \, {\rm M}_{\odot}$, mergers have the expected effect: if gas-poor they promote the formation of spheroidal galaxies, whereas gas-rich mergers favour the formation and survivability of massive discs. This trend, however, breaks at lower masses. For objects with $M_{\ast} \lesssim 10^{11} \, {\rm M}_{\odot}$, mergers do not seem to play any significant role in determining the morphology, with accreted stellar fractions and mean merger gas fractions that are indistinguishable between spheroidal and disc-dominated galaxies. On the other hand, halo spin correlates with morphology primarily in the least massive objects in the sample ($M_{\ast} \lesssim 10^{10} \, {\rm M}_{\odot}$), but only weakly for galaxies above that mass. Our results support a scenario where (1) mergers play a dominant role in shaping the morphology of massive galaxies, (2) halo spin is important for the morphology of dwarfs, and (3) the morphology of medium-sized galaxies -- including the Milky Way -- shows little dependence on galaxy assembly history or halo spin, at least when these two factors are considered individually.

We use the Illustris simulation to study the relative contributions of in situ star formation and stellar accretion to the build-up of galaxies over an unprecedentedly wide range of masses ($M_{\ast} = 10^9-10^{12} \, {\rm M_{\odot}}$), galaxy types, environments, and assembly histories. We find that the `two-phase' picture of galaxy formation predicted by some models is a good approximation only for the most massive galaxies in our simulation -- namely, the stellar mass growth of galaxies below a few times $10^{11} \, {\rm M_{\odot}}$ is dominated by in situ star formation at all redshifts. The fraction of the total stellar mass of galaxies at $z=0$ contributed by accreted stars shows a strong dependence on galaxy stellar mass, ranging from about 10 per cent for Milky Way-sized galaxies to over 80 per cent for $M_{\ast} \approx 10^{12} \, {\rm M_{\odot}}$ objects, yet with a large galaxy-to-galaxy variation. At a fixed stellar mass, elliptical galaxies and those formed at the centres of younger haloes exhibit larger fractions of ex situ stars than disc-like galaxies and those formed in older haloes. On average, $\sim$50 per cent of the ex situ stellar mass comes from major mergers (stellar mass ratio $\mu > 1/4$), $\sim$20 per cent from minor mergers ($1/10 < \mu < 1/4$), $\sim$20 per cent from very minor mergers ($\mu < 1/10$), and $\sim$10 per cent from stars that were stripped from surviving galaxies (e.g. flybys or ongoing mergers). These components are spatially segregated, with in situ stars dominating the innermost regions of galaxies, and ex situ stars being deposited at larger galactocentric distances in order of decreasing merger mass ratio.

Feedback from central supermassive blackholes is often invoked to explain the low star formation rates in massive galaxies at the centers of galaxy clusters. However, the detailed physics of the coupling of the injected feedback energy with the intracluster medium is still unclear. Using high-resolution magnetohydrodynamic cosmological simulations of galaxy cluster formation, we investigate the role of anisotropic thermal conduction in shaping the thermodynamic structure of clusters, and, in particular, in modifying the impact of black hole feedback. Stratified anisotropically conducting plasmas are formally always unstable, and thus more prone to mixing, an expectation borne out by our results. The increased mixing efficiently isotropizes the injected feedback energy which in turn significantly improves the coupling between the feedback energy and the intracluster medium. This facilitates an earlier disruption of the cool core, reduces the star formation rate by more than an order of magnitude, and results in earlier quenching despite an overall lower amount of feedback energy injected into the cluster core. With conduction, the metallicity gradients and dispersions are lowered, aligning them better with observational constraints. These results highlight the important role of thermal conduction in establishing and maintaining quiescence of massive galaxies.

We investigate the effects of self-interacting dark matter (SIDM) on the tidal stripping and evaporation of satellite galaxies in a Milky Way-like host. We use a suite of five zoom-in, dark-matter-only simulations, two with velocity-independent SIDM cross sections, two with velocity-dependent SIDM cross sections, and one cold dark matter simulation for comparison. After carefully assigning stellar mass to satellites at infall, we find that stars are stripped at a higher rate in SIDM than in CDM. In contrast, the total bound dark matter mass loss rate is minimally affected, with subhalo evaporation having negligible effects on satellites for viable SIDM models. Centrally located stars in SIDM haloes disperse out to larger radii as cores grow. Consequently, the half-light radius of satellites increases, stars become more vulnerable to tidal stripping, and the stellar mass function is suppressed. We find that the ratio of core radius to tidal radius accurately predicts the relative strength of enhanced SIDM stellar stripping. Velocity-independent SIDM models show a modest increase in the stellar stripping effect with satellite mass, whereas velocity-dependent SIDM models show a large increase in this effect towards lower masses, making observations of ultra-faint dwarfs prime targets for distinguishing between and constraining SIDM models. Due to small cores in the largest satellites of velocity-dependent SIDM, no identifiable imprint is left on the all-sky properties of the stellar halo. While our results focus on SIDM, the main physical mechanism of enhanced tidal stripping of stars apply similarly to satellites with cores formed via other means.

We present a constrained transport (CT) algorithm for solving the 3D ideal magnetohydrodynamic (MHD) equations on a moving mesh, which maintains the divergence-free condition on the magnetic field to machine-precision. Our CT scheme uses an unstructured representation of the magnetic vector potential, making the numerical method simple and computationally efficient. The scheme is implemented in the moving mesh code Arepo. We demonstrate the performance of the approach with simulations of driven MHD turbulence, a magnetized disc galaxy, and a cosmological volume with primordial magnetic field. We compare the outcomes of these experiments to those obtained with a previously implemented Powell divergence-cleaning scheme. While CT and the Powell technique yield similar results in idealized test problems, some differences are seen in situations more representative of astrophysical flows. In the turbulence simulations, the Powell cleaning scheme artificially grows the mean magnetic field, while CT maintains this conserved quantity of ideal MHD. In the disc simulation, CT gives slower magnetic field growth rate and saturates to equipartition between the turbulent kinetic energy and magnetic energy, whereas Powell cleaning produces a dynamically dominant magnetic field. Such difference has been observed in adaptive-mesh refinement codes with CT and smoothed-particle hydrodynamics codes with divergence-cleaning. In the cosmological simulation, both approaches give similar magnetic amplification, but Powell exhibits more cell-level noise. CT methods in general are more accurate than divergence-cleaning techniques, and, when coupled to a moving mesh can exploit the advantages of automatic spatial/temporal adaptivity and reduced advection errors, allowing for improved astrophysical MHD simulations.

The velocity function is a fundamental observable statistic of the galaxy population, similarly impor- tant as the luminosity function, but much more difficult to measure. In this work we present the first directly measured circular velocity function that is representative between 60 < v_circ < 320 km/s for galaxies of all morphological types at a given rotation velocity. For the low mass galaxy population (60 < v_circ < 170 km/s), we use the HIPASS velocity function. For the massive galaxy population (170 < v_circ < 320 km/s), we use stellar circular velocities from the Calar Alto Legacy Integral Field Area Survey (CALIFA). In earlier work we obtained the measurements of circular velocity at the 80% light radius for 226 galaxies and demonstrated that the CALIFA sample can produce volume- corrected galaxy distribution functions. The CALIFA velocity function includes homogeneous velocity measurements of both late and early-type rotation-supported galaxies and has the crucial advantage of not missing gas-poor massive ellipticals that HI surveys are blind to. We show that both velocity functions can be combined in a seamless manner, as their ranges of validity overlap. The resulting observed velocity function is compared to velocity functions derived from cosmological simulations of the z = 0 galaxy population. We find that dark matter-only simulations show a strong mismatch with the observed VF. Hydrodynamic simulations fare better, but still do not fully reproduce observations.

We use particle data from the Illustris simulation, combined with individual kinematic constraints on the mass of the Milky Way (MW) at specific distances from the Galactic center, to infer the radial distribution of the MW's dark matter halo mass. Our method allows us to convert any constraint on the mass of the MW within a fixed distance to a full circular velocity profile to the MW's virial radius. As primary examples, we take two recent (and discrepant) measurements of the total mass within 50 kpc of the Galaxy and find they imply very different mass profiles and stellar masses for the Galaxy. The dark-matter-only version of the Illustris simulation enables us to compute the effects of galaxy formation on such constraints on a halo-by-halo basis; on small scales, galaxy formation enhances the density relative to dark-matter-only runs, while the total mass density is approximately 20% lower at large Galactocentric distances. We are also able to quantify how current and future constraints on the mass of the MW within specific radii will be reflected in uncertainties on its virial mass: even a measurement of M(<50 kpc) with essentially perfect precision still results in a 20% uncertainty on the virial mass of the Galaxy, while a future measurement of M(<100 kpc) with 10% errors would result in the same level of uncertainty. We expect that our technique will become even more useful as (1) better kinematic constraints become available at larger distances and (2) cosmological simulations provide even more faithful representations of the observable Universe.

Hydrodynamical shocks are a manifestation of the non-linearity of the Euler equations and play a fundamental role in cosmological gas dynamics. In this work, we identify and analyse shocks in the Illustris simulation, and contrast the results with those of non-radiative runs. We show that simulations with more comprehensive physical models of galaxy formation pose new challenges for shock finding algorithms due to radiative cooling and star-forming processes, prompting us to develop a number of methodology improvements. We find in Illustris a total shock surface area which is about 1.4 times larger at the present epoch compared to non-radiative runs, and an energy dissipation rate at shocks which is higher by a factor of around 7. Remarkably, shocks with Mach numbers above and below $\mathcal{M}\approx10$ contribute about equally to the total dissipation across cosmic time. This is in sharp contrast to non-radiative simulations, and we demonstrate that a large part of the difference arises due to strong black hole radio-mode feedback in Illustris. We also provide an overview of the large diversity of shock morphologies, which includes complex networks of halo-internal shocks, shocks on to cosmic sheets, feedback shocks due to black holes and galactic winds, as well as ubiquitous accretion shocks. In high redshift systems more massive than $10^{12}\,\mathrm{M}_\odot$ we discover the existence of a double accretion shock pattern in haloes. They are created when gas streams along filaments without being shocked at the outer accretion shock, but then forms a second, roughly spherical accretion shock further inside.

Galaxy comoving number density is commonly used to forge progenitor/descendant links between observed galaxy populations at different epochs. However, this method breaks down in the presence of galaxy mergers, or when galaxies experience stochastic growth rates. We present a simple analytic framework to treat the physical processes that drive the evolution and diffusion of galaxies within comoving number density space. The evolution in mass rank order of a galaxy population with time is influenced by the galaxy coagulation rate and galaxy "mass rank scatter" rate. We quantify the relative contribution of these two effects to the mass rank order evolution. We show that galaxy coagulation is dominant at lower redshifts and stellar masses, while scattered growth rates dominate the mass rank evolution at higher redshifts and stellar masses. For a galaxy population at $10^{10} M_\odot$, coagulation has been the dominant effect since $z=2.2$, but a galaxy population at $10^{11} M_\odot$ was dominated by mass rank scatter until $z=0.6$. We show that although the forward and backward median number density evolution tracks are asymmetric, the backward median number density evolution can be obtained by convolving the descendant distribution function with progenitor relative abundances. We tabulate fits for the median number density evolution and scatter which can be applied to improve the way galaxy populations are linked in multi-epoch observational datasets.

We present the first simulations within an effective theory of structure formation (ETHOS), which includes the effect of interactions between dark matter and dark radiation on the linear initial power spectrum and dark matter self-interactions during non-linear structure formation. We simulate a Milky Way-like halo in four different dark matter models and the cold dark matter case. Our highest resolution simulation has a particle mass of $2.8\times 10^4\,{\rm M}_\odot$ and a softening length of $72.4\,{\rm pc}$. We demonstrate that all alternative models have only a negligible impact on large scale structure formation. On galactic scales, however, the models significantly affect the structure and abundance of subhaloes due to the combined effects of small scale primordial damping in the power spectrum and late time self-interactions. We derive an analytic mapping from the primordial damping scale in the power spectrum to the cutoff scale in the halo mass function and the kinetic decoupling temperature. We demonstrate that certain models within this extended effective framework that can alleviate the too-big-to-fail and missing satellite problems simultaneously, and possibly the core-cusp problem. The primordial power spectrum cutoff of our models naturally creates a diversity in the circular velocity profiles, which is larger than that found for cold dark matter simulations. We show that the parameter space of models can be constrained by contrasting model predictions to astrophysical observations. For example, some models may be challenged by the missing satellite problem if baryonic processes were to be included and even over-solve the too-big-to-fail problem; thus ruling them out.

The concept of the cosmic web, viewing the Universe as a set of discrete galaxies held together by gravity, is deeply engrained in cosmology. Yet, little is known about the most effective construction and the characteristics of the underlying network. Here we explore seven network construction algorithms that use various galaxy properties, from their location, to their size and relative velocity, to assign a network to galaxy distributions provided by both simulations and observations. We find that a model relying only on spatial proximity offers the best correlations between the physical characteristics of the connected galaxies. We show that the properties of the networks generated from simulations and observations are identical, unveiling a deep universality of the cosmic web.

(abridged) The large-scale distribution of galaxies is generally analyzed using the two-point correlation function. However, this statistic does not capture the topology of the distribution, and it is necessary to resort to higher order correlations to break degeneracies. We demonstrate that an alternate approach using network analysis can discriminate between topologically different distributions that have similar two-point correlations. We investigate two galaxy point distributions, one produced by a cosmological simulation and the other by a L\'evy walk. For the cosmological simulation, we adopt the redshift $z = 0.58$ slice from Illustris (Vogelsberger et al. 2014A) and select galaxies with stellar masses greater than $10^8$$M_\odot$. The two point correlation function of these simulated galaxies follows a single power-law, $\xi(r) \sim r^{-1.5}$. Then, we generate L\'evy walks matching the correlation function and abundance with the simulated galaxies. We find that, while the two simulated galaxy point distributions have the same abundance and two point correlation function, their spatial distributions are very different; most prominently, \emph{filamentary structures}, absent in L\'evy fractals. To quantify these missing topologies, we adopt network analysis tools and measure diameter, giant component, and transitivity from networks built by a conventional friends-of-friends recipe with various linking lengths. Unlike the abundance and two point correlation function, these network quantities reveal a clear separation between the two simulated distributions; therefore, the galaxy distribution simulated by Illustris is not a L\'evy fractal quantitatively. We find that the described network quantities offer an efficient tool for discriminating topologies and for comparing observed and theoretical distributions.

It is widely believed that the large discrepancy between the observed number of satellite galaxies and the predicted number of dark subhalos can be resolved via a variety of baryonic effects which suppress star formation in low mass halos.Supporting this hypothesis, numerous high resolution simulations with star formation, and associated feedback have been shown to reproduce the satellite luminosity function around Milky Way-mass simulated galaxies at redshift zero. However, a more stringent test of these models is their ability to simultaneously match the satellite luminosity functions of a range of host halo masses and redshifts. In this work we measure the luminosity function of faint (sub-Small Magellanic Cloud luminosity) satellites around hosts with stellar masses 10.5$<\log_{10}$M$_*$/M$_\odot<11.5$ to an unprecedented redshift of 1.5. This new measurement of the satellite luminosity function provides powerful new constraining power; we compare these results with predictions from four different simulations and show that although the models perform similarly over-all, no one model reproduces the satellite luminosity function reliably at all redshifts and host stellar masses. This result highlights the continued need for improvement in understanding the fundamental physics that governs satellite galaxy evolution.

We introduce a dust model for cosmological simulations implemented in the moving-mesh code AREPO and present a suite of cosmological hydrodynamical zoom-in simulations to study dust formation within galactic haloes. Our model accounts for the stellar production of dust, accretion of gas-phase metals onto existing grains, destruction of dust through local supernova activity, and dust driven by winds from star-forming regions. We find that accurate stellar and active galactic nuclei feedback is needed to reproduce the observed dust-metallicity relation and that dust growth largely dominates dust destruction. Our simulations predict a dust content of the interstellar medium which is consistent with observed scaling relations at $z = 0$, including scalings between dust-to-gas ratio and metallicity, dust mass and gas mass, dust-to-gas ratio and stellar mass, and dust-to-stellar mass ratio and gas fraction. We find that roughly two-thirds of dust at $z = 0$ originated from Type II supernovae, with the contribution from asymptotic giant branch stars below 20 per cent for $z \gtrsim 5$. While our suite of Milky Way-sized galaxies forms dust in good agreement with a number of key observables, it predicts a high dust-to-metal ratio in the circumgalactic medium, which motivates a more realistic treatment of thermal sputtering of grains and dust cooling channels.

Observations at low redshifts thus far fail to account for all of the baryons expected in the Universe according to cosmological constraints. A large fraction of the baryons presumably resides in a thin and warm-hot medium between the galaxies, where they are difficult to observe due to their low densities and high temperatures. Cosmological simulations of structure formation can be used to verify this picture and provide quantitative predictions for the distribution of mass in different large-scale structure components. Here we study the distribution of baryons and dark matter at different epochs using data from the Illustris simulation. We identify regions of different dark matter density with the primary constituents of large-scale structure, allowing us to measure mass and volume of haloes, filaments and voids. At redshift zero, we find that 49 % of the dark matter and 23 % of the baryons are within haloes more massive than the resolution limit of $2\times 10^8$ M$_\odot$. The filaments of the cosmic web host a further 45 % of the dark matter and 46 % of the baryons. The remaining 31 % of the baryons reside in voids. The majority of these baryons have been transported there through active galactic nuclei feedback. We note that the feedback model of Illustris is too strong for heavy haloes, therefore it is likely that we are overestimating this amount. Categorizing the baryons according to their density and temperature, we find that 17.8 % of them are in a condensed state, 21.6 % are present as cold, diffuse gas, and 53.9 % are found in the state of a warm-hot intergalactic medium.

We present the largest number of Milky Way sized dark matter halos simulated at very high mass ($\sim$$10^4$ M$_\odot$/particle) and temporal resolution ($\sim$5 Myrs/snapshot) done to date, quadrupling what is currently available in the literature. This initial suite consists of the first 24 halos of the $Caterpillar$ $Project$ (www.caterpillarproject.org) whose project goal of 60 - 70 halos will be made public when complete. We resolve $\sim$20,000 gravitationally bound subhalos within the virial radius of each host halo. Over the ranges set by our spatial resolution our convergence is excellent and improvements were made upon current state-of-the-art halo finders to better identify substructure at such high resolutions (e.g., on average we recover $\sim$4 subhalos in each host halo above 10$^8$ M$_\odot$ which would have otherwise not been found using conventional methods). For our relaxed halos, the inner profiles are reasonably fit by Einasto profiles ($\alpha$ = 0.169 $\pm$ 0.023) though this depends on the relaxed nature and assembly history of a given halo. Averaging over all halos, the substructure mass fraction is $f_{m,subs} = 0.121 \pm 0.041$, and mass function slope is d$N$/d$M\propto M^{-1.88 \pm 0.10}$ though we find scatter in the normalizations for fixed halo mass due to more concentrated hosts having less subhalos at fixed subhalo mass. There are no biases stemming from Lagrangian volume selection as all Lagrangian volume types are included in our sample. Our detailed contamination study of 264 low resolution halos has resulted in obtaining very large and unprecedented, high-resolution regions around our host halos for our target resolution (sphere of radius $\sim$$1.4 \pm 0.4$ Mpc) allowing for accurate studies of low mass dwarf galaxies at large galactocentric radii and the very first stellar systems at high redshift ($z \geq$ 10).

We examine the stability of feedback-regulated star formation (SF) in galactic nuclei and contrast it to SF in extended discs. In galactic nuclei the dynamical time becomes shorter than the time over which feedback from young stars evolves. We argue analytically that the balance between stellar feedback and gravity is unstable in this regime. We study this using numerical simulations with pc-scale resolution and explicit stellar feedback taken from stellar evolution models. The nuclear gas mass, young stellar mass, and SFR within the central ~100 pc (the short-timescale regime) never reach steady-state, but instead go through dramatic, oscillatory cycles. Stars form until a critical surface density of young stars is present (such that feedback overwhelms gravity), at which point they begin to expel gas from the nucleus. Since the dynamical times are shorter than the stellar evolution times, the stars do not die as the gas is expelled, but continue to push, triggering a runaway quenching of star formation in the nucleus. However the expelled gas is largely not unbound from the galaxy, but goes into a galactic fountain which re-fills the nuclear region after the massive stars from the previous burst cycle have died off (~50 Myr timescale). On large scales (>1 kpc), the galaxy-scale gas content and SFR is more stable. We examine the consequences of this episodic nuclear star formation for the Kennicutt-Schmidt (KS) relation: while a tight KS relation exists on ~1 kpc scales in good agreement with observations, the scatter increases dramatically in smaller apertures centered on galactic nuclei.

Using adiabatic high-resolution numerical simulations we quantify the effect of the streaming motion of baryons with respect to dark matter at the time of recombination on structure formation and evolution. Formally a second order effect, the baryonic stream velocity has proven to have significant impact on dark matter halo abundance, as well as on the gas content and morphology of small galaxy clusters. In this work, we study the impact of stream velocity on the formation and gas content of haloes with masses up to $10^9 M_{\odot}$, an order of magnitude larger than previous studies. We find that the non-zero stream velocity has a sizable impact on the number density of haloes with masses $\lesssim$ few $\times 10^7 M_{\odot}$ up to $z=10$, the final redshift of our simulations. Furthermore, the gas stream velocity induces a suppression of the gas fraction in haloes, which at z=10 is $\sim 10 \%$ for objects with $M\sim10^7M_{\odot}$, as well as a flattening of the gas density profiles in the inner regions of haloes. We further identify and study the formation, in the context of a non-zero stream velocity, of moderately long lived gas dominated structures at intermediate redshifts $10 < z < 20$, which Naoz and Narayan have recently proposed as potential progenitors of globular clusters.

Supermassive black hole (BH) mergers produce powerful gravitational wave (GW) emission. Asymmetry in this emission imparts a recoil kick to the merged BH, which can eject the BH from its host galaxy altogether. Recoiling BHs could be observed as offset active galactic nuclei (AGN). Several candidates have been identified, but systematic searches have been hampered by large uncertainties regarding their observability. By extracting merging BHs and host galaxy properties from the Illustris cosmological simulations, we have developed a comprehensive model for recoiling AGN. Here, for the first time, we model the effects of BH spin alignment and recoil dynamics based on the gas-richness of host galaxies. We predict that if BH spins are not highly aligned, seeing-limited observations could resolve offset AGN, making them promising targets for all-sky surveys. For randomly-oriented spins, less than about 10 spatially-offset AGN may be detectable in HST-COSMOS, and > 10^3 could be found with Pan-STARRS, LSST, Euclid, and WFIRST. Nearly a thousand velocity-offset AGN are predicted within the SDSS footprint; the rarity of large broad-line offsets among SDSS quasars is likely due in part to selection effects but suggests that spin alignment plays a role in suppressing recoils. Nonetheless, in our most physically motivated model where alignment occurs only in gas-rich mergers, hundreds of offset AGN should be found in all-sky surveys. Our findings strongly motivate a dedicated search for recoiling AGN.

We present a novel implementation of an extremum preserving anisotropic diffusion solver for thermal conduction on the unstructured moving Voronoi mesh of the AREPO code. The method relies on splitting the one-sided facet fluxes into normal and oblique components, with the oblique fluxes being limited such that the total flux is both locally conservative and extremum preserving. The approach makes use of harmonic averaging points and a simple, robust interpolation scheme that works well for strong heterogeneous and anisotropic diffusion problems. Moreover, the required discretisation stencil is small. Efficient fully implicit and semi-implicit time integration schemes are also implemented. We perform several numerical tests that evaluate the stability and accuracy of the scheme, including applications such as point explosions with heat conduction and calculations of convective instabilities in conducting plasmas. The new implementation is suitable for studying important astrophysical phenomena, such as the conductive heat transport in galaxy clusters, the evolution of supernova remnants, or the distribution of heat from blackhole-driven jets into the intracluster medium.

Massive quiescent galaxies have much smaller physical sizes at high redshift than today. The strong evolution of galaxy size may be caused by progenitor bias, major and minor mergers, adiabatic expansion, and/or renewed star formation, but it is difficult to test these theories observationally. Herein, we select a sample of 35 massive, compact galaxies ($M_* = 1-3 \times 10^{11}$ M$_\odot$, $M_*/R^{1.5} > 10^{10.5}$ M$_\odot$/kpc$^{1.5}$) at $z=2$ in the cosmological hydrodynamical simulation Illustris and trace them forward to $z=0$ to uncover their evolution and identify their descendants. By $z=0$, the original factor of 3 difference in stellar mass spreads to a factor of 20. The dark matter halo masses similarly spread from a factor of 5 to 40. The galaxies' evolutionary paths are diverse: about half acquire an ex-situ envelope and are the core of a more massive descendant, a third survive undisturbed and gain very little mass, 15% are consumed in a merger with a more massive galaxy, and a small remainder are thoroughly mixed by major mergers. The galaxies grow in size as well as mass, and only $\sim$10% remain compact by $z=0$. The majority of the size growth is driven by the acquisition of ex-situ mass. The most massive galaxies at $z=0$ are the most likely to have compact progenitors, but this trend possesses significant dispersion which precludes a direct linkage to compact galaxies at $z=2$. The compact galaxies' merger rates are influenced by their $z=2$ environments, so that isolated or satellite compact galaxies (which are protected from mergers) are the most likely to survive to the present day.

We investigate the effects of varying the intensity of the primordial magnetic seed field on the global properties of the galaxy population in ideal magnetohydrodynamic cosmological simulations performed with the moving-mesh code AREPO. We vary the seed field in our calculations in a range of values still compatible with the current cosmological upper limits. We show that above a critical intensity of $\simeq 10^{-9}\,{\rm G}$, the additional pressure arising from the field strongly affects the evolution of gaseous structures, leading to a suppression of the cosmic star formation history, \rev{which is stronger for larger seed fields. This directly reflects into a lower total galaxy count above a fixed stellar mass threshold at all redshifts, and} a lower galaxy number density at fixed stellar mass and a less massive stellar component at fixed virial mass at all mass scales. These signatures may be used, in addition to the existing methods, to derive tighter constraints on primordial magnetic seed field intensities.

Galaxy clusters contain a large population of low mass dwarf elliptical galaxies whose exact origin is unclear: their colors, structural properties and kinematics differ substantially from those of dwarf irregulars in the field. We use the Illustris cosmological simulation to study differences in the assembly paths of dwarf galaxies (3e8 < M_*/M_sun < 1e10) according to their environment. We find that cluster dwarfs achieve their maximum total and stellar mass on average ~ 8 and ~ 4.5 Gyr ago (or redshifts z = 1.0 and z = 0.4, respectively), around the time of infall into the clusters. In contrast, field dwarfs not subjected to environmental stripping, reach their maximum mass at redshift z = 0. This different assembly history naturally produces a color bimodality, with blue isolated dwarfs and redder cluster dwarfs exhibiting negligible star-formation today. The cessation of star formation happens over median times 3.5-5 Gyr depending on stellar mass, and shows a large scatter (~ 1-8 Gyr), with the lower values associated with starburst events that occur at infall through the virial radius or pericentric passages. We argue that such starbursts together with the early assembly of cluster dwarfs can provide a natural explanation for the higher specific frequency of globular clusters (GCs) in cluster dwarfs, as found observationally. We present a simple model for the formation and stripping of GCs that supports this interpretation. The origin of dwarf ellipticals in clusters is, therefore, consistent with an environmentally-driven evolution of field dwarf irregulars. However the z = 0 field analogs of cluster dwarf progenitors have today stellar masses a factor ~ 3 larger --a difference arising from the early truncation of star formation in cluster dwarfs.

We analyze the mass, temperature, metal enrichment, and OVI abundance of the circumgalactic medium (CGM) around $z\sim 0.2$ galaxies of mass $10^9 M_\odot <M_\bigstar < 10^{11.5} M_\odot$ in the Illustris simulation. Among star-forming galaxies, the mass, temperature, and metallicity of the CGM increase with stellar mass, driving an increase in the OVI column density profile of $\sim 0.5$ dex with each $0.5$ dex increase in stellar mass. Observed OVI column density profiles exhibit a weaker mass dependence than predicted: the simulated OVI abundance profiles are consistent with those observed for star-forming galaxies of mass $M_\bigstar = 10^{10.5-11.5} M_\odot$, but underpredict the observed OVI abundances by $\gtrsim 0.8$ dex for lower-mass galaxies. We suggest that this discrepancy may be alleviated with additional heating of the abundant cool gas in low-mass halos, or with increased numerical resolution capturing turbulent/conductive mixing layers between CGM phases. Quenched galaxies of mass $M_\bigstar = 10^{10.5-11.5} M_\odot$ are found to have 0.3-0.8 dex lower OVI column density profiles than star-forming galaxies of the same mass, in qualitative agreement with the observed OVI abundance bimodality. This offset is driven by AGN feedback, which quenches galaxies by heating the CGM and ejecting significant amounts of gas from the halo. Finally, we find that the inclusion of the central galaxy's radiation field may enhance the photoionization of the CGM within $\sim 50$ kpc, further increasing the predicted OVI abundance around star-forming galaxies.

We present the full public release of all data from the Illustris simulation project. Illustris is a suite of large volume, cosmological hydrodynamical simulations run with the moving-mesh code Arepo and including a comprehensive set of physical models critical for following the formation and evolution of galaxies across cosmic time. Each simulates a volume of (106.5 Mpc)^3 and self-consistently evolves five different types of resolution elements from a starting redshift of z=127 to the present day, z=0. These components are: dark matter particles, gas cells, passive gas tracers, stars and stellar wind particles, and supermassive black holes. This data release includes the snapshots at all 136 available redshifts, halo and subhalo catalogs at each snapshot, and two distinct merger trees. Six primary realizations of the Illustris volume are released, including the flagship Illustris-1 run. These include three resolution levels with the fiducial "full" baryonic physics model, and a dark matter only analog for each. In addition, we provide four distinct, high time resolution, smaller volume "subboxes". The total data volume is ~265 TB, including ~800 full volume snapshots and ~30,000 subbox snapshots. We describe the released data products as well as tools we have developed for their analysis. All data may be directly downloaded in its native HDF5 format. Additionally, we release a comprehensive, web-based API which allows programmatic access to search and data processing tasks. In both cases we provide example scripts and a getting-started guide in several languages: currently, IDL, Python, and Matlab. This paper addresses scientific issues relevant for the interpretation of the simulations, serves as a pointer to published and on-line documentation of the project, describes planned future additional data releases, and discusses technical aspects of the release.

We study how optical galaxy morphology depends on mass and star formation rate (SFR) in the Illustris Simulation. To do so, we measure automated galaxy structures in 10808 simulated galaxies at z=0 with stellar masses 10^9.7 < M_*/M_sun < 10^12.3. We add observational realism to idealized synthetic images and measure non-parametric statistics in rest-frame optical and near-IR images from four directions. We find that Illustris creates a morphologically diverse galaxy population, occupying the observed bulge strength locus and reproducing median morphology trends versus stellar mass, SFR, and compactness. Morphology correlates realistically with rotation, following classification schemes put forth by kinematic surveys. Type fractions as a function of environment agree roughly with data. These results imply that connections among mass, star formation, and galaxy structure arise naturally from models matching global star formation and halo occupation functions when simulated with accurate methods. This raises a question of how to construct experiments on galaxy surveys to better distinguish between models. We predict that at fixed halo mass near 10^12 M_sun, disc-dominated galaxies have higher stellar mass than bulge-dominated ones, a possible consequence of the Illustris feedback model. While Illustris galaxies at M_* ~ 10^11 M_sun have a reasonable size distribution, those at M_* ~ 10^10 M_sun have half-light radii larger than observed by a factor of two. Furthermore, at M_* ~ 10^10.5-10^11 M_sun, a relevant fraction of Illustris galaxies have distinct "ring-like" features, such that the bright pixels have an unusually wide spatial extent.

Comparisons between observational surveys and galaxy formation models find that the mass of dark matter haloes can largely explain galaxies' stellar mass. However, it remains uncertain whether additional environmental variables, generally referred to as assembly bias, are necessary to explain other galaxy properties. We use the Illustris Simulation to investigate the role of assembly bias in producing galactic conformity by considering 18,000 galaxies with $M_{stellar}$ > $2 \times 10^9$ $M_{\odot}$. We find a significant signal of galactic conformity: out to distances of about 10 Mpc, the mean red fraction of galaxies around redder galaxies is higher than around bluer galaxies at fixed stellar mass. Dark matter haloes exhibit an analogous conformity signal, in which the fraction of haloes formed at earlier times (old haloes) is higher around old haloes than around younger ones at fixed halo mass. A plausible interpretation of galactic conformity can be given as a combination of the halo conformity signal with the galaxy color-halo age relation: at fixed stellar mass, particularly toward the low-mass end, Illustris' galaxy colors correlate with halo age, with the reddest galaxies (often satellites) being preferentially found in the oldest haloes. In fact, we can explain the galactic conformity effect with a simple semi-empirical model, by assigning stellar mass based on halo mass (abundance matching) and by assigning galaxy color based on halo age (age matching). We investigate other interpretations for the galactic conformity, particularly its dependence on the isolation criterion and on the central-satellite information. Regarding comparison to observations, we conclude that the adopted selection/isolation criteria, projection effects, and stacking techniques can have a significant impact on the measured amplitude of the conformity signal.

We perform uniformly sampled large-scale cosmological simulations including magnetic fields with the moving mesh code AREPO. We run two sets of MHD simulations: one including adiabatic gas physics only; the other featuring the fiducial feedback model of the Illustris simulation. In the adiabatic case, the magnetic field amplification follows the $B \propto \rho^{2/3}$ scaling derived from `flux-freezing' arguments, with the seed field strength providing an overall normalization factor. At high baryon overdensities the amplification is enhanced by shear flows and turbulence. Feedback physics and the inclusion of radiative cooling change this picture dramatically. In haloes, gas collapses to much larger densities and the magnetic field is amplified strongly and to the same maximum intensity irrespective of the initial seed field of which any memory is lost. At lower densities a dependence on the seed field strength and orientation, which in principle can be used to constrain models of cosmic magnetogenesis, is still present. Inside the most massive haloes magnetic fields reach values of $\sim 10-100\,\,{\rm \mu G}$, in agreement with galaxy cluster observations. The topology of the field is tangled and gives rise to rotation measure signals in reasonable agreement with the observations. However, the rotation measure signal declines too rapidly towards larger radii as compared to observational data.

Hydrodynamical simulations of galaxy formation such as the Illustris simulations have progressed to a state where they approximately reproduce the observed stellar mass function from high to low redshift. This in principle allows self-consistent models of reionization that exploit the accurate representation of the diffuse gas distribution together with the realistic growth of galaxies provided by these simulations, within a representative cosmological volume. In this work, we apply and compare two radiative transfer algorithms implemented in a GPU-accelerated code to the $106.5\,{\rm Mpc}$ wide volume of Illustris in postprocessing in order to investigate the reionization transition predicted by this model. We find that the first generation of galaxies formed by Illustris is just about able to reionize the universe by redshift $z\sim 7$, provided quite optimistic assumptions about the escape fraction and the resolution limitations are made. Our most optimistic model finds an optical depth of $\tau\simeq 0.065$, which is in very good agreement with recent Planck 2015 determinations. Furthermore, we show that moment-based approaches for radiative transfer with the M1 closure give broadly consistent results with our angular-resolved radiative transfer scheme. In our favoured fiducial model, 20\% of the hydrogen is reionized by redshift $z=9.20$, and this rapidly climbs to 80\% by redshift $z=6.92$. It then takes until $z=6.24$ before 99\% of the hydrogen is ionized. On average, reionization proceeds `inside-out' in our models, with a size distribution of reionized bubbles that progressively features regions of ever larger size while the abundance of small bubbles stays fairly constant.

The cumulative comoving number-density of galaxies as a function of stellar mass or central velocity dispersion is commonly used to link galaxy populations across different epochs. By assuming that galaxies preserve their number-density in time, one can infer the evolution of their properties, such as masses, sizes, and morphologies. However, this assumption does not hold in the presence of galaxy mergers or when rank ordering is broken owing to variable stellar growth rates. We present an analysis of the evolving comoving number density of galaxy populations found in the Illustris cosmological hydrodynamical simulation focused on the redshift range $0\leq z \leq 3$. Our primary results are as follows: 1) The inferred average stellar mass evolution obtained via a constant comoving number density assumption is systematically biased compared to the merger tree results at the factor of $\sim$2(4) level when tracking galaxies from redshift $z=0$ out to redshift $z=2(3)$; 2) The median number density evolution for galaxy populations tracked forward in time is shallower than for galaxy populations tracked backward in time; 3) A similar evolution in the median number density of tracked galaxy populations is found regardless of whether number density is assigned via stellar mass, stellar velocity dispersion, or dark matter halo mass; 4) Explicit tracking reveals a large diversity in galaxies' assembly histories that cannot be captured by constant number-density analyses; 5) The significant scatter in galaxy linking methods is only marginally reduced by considering a number of additional physical and observable galaxy properties as realized in our simulation. We provide fits for the forward and backward median evolution in stellar mass and number density and discuss implications of our analysis for interpreting multi-epoch galaxy property observations.

We study the properties of black holes and their host galaxies across cosmic time in the Illustris simulation. Illustris is a large scale cosmological hydrodynamical simulation which resolves a (106.5 Mpc)^3 volume with more than 12 billion resolution elements and includes state-of-the-art physical models relevant for galaxy formation. We find that the black hole mass density for redshifts z = 0 - 5 and the black hole mass function at z = 0 predicted by Illustris are in very good agreement with the most recent observational constraints. We show that the bolometric and hard X-ray luminosity functions of AGN at z = 0 and 1 reproduce observational data very well over the full dynamic range probed. Unless the bolometric corrections are largely underestimated, this requires radiative efficiencies to be on average low, epsilon_r <= 0.1, noting however that in our model radiative efficiencies are degenerate with black hole feedback efficiencies. Cosmic downsizing of the AGN population is in broad agreement with the findings from X-ray surveys, but we predict a larger number density of faint AGN at high redshifts than currently inferred. We also study black hole -- host galaxy scaling relations as a function of galaxy morphology, colour and specific star formation rate. We find that black holes and galaxies co-evolve at the massive end, but for low mass, blue and star-forming galaxies there is no tight relation with either their central black hole masses or the nuclear AGN activity.

We study how X-rays from stellar binary systems and the hot intracluster medium (ICM) affect the radiative cooling rates of gas in galaxies. Our study uses a novel implementation of gas cooling in the moving-mesh hydrodynamics code \textsc{arepo}. X-rays from stellar binaries do not affect cooling at all as their emission spectrum is too hard to effectively couple with galactic gas. In contrast, X-rays from the ICM couple well with gas in the temperature range $10^4 - 10^6$ K. Idealised simulations show that the hot halo radiation field has minimal impact on the dynamics of cooling flows in clusters because of the high virial temperature ($> 10^7$K), making the interaction between the gas and incident photons very ineffective. Satellite galaxies in cluster environments, on the other hand, experience a high radiation flux due to the emission from the host halo. Low mass satellites ($< 10^{12}\rm{M_\odot}$) in particular have virial temperatures that are exactly in the regime where the effect of the radiation field is maximal. Idealised simulations of satellite galaxies including only the effect of host halo radiation (no ram pressure stripping or tidal effects) fields show a drastic reduction in the amount of cool gas formed ($\sim 40\%$) on a short timescale of about $0.5$ Gyrs. A galaxy merger simulation including all the other environmental quenching mechanisms, shows about $20\%$ reduction in the stellar mass of the satellite and about $\sim 30\%$ reduction in star formation rate after $1$ Gyr due to the host hot halo radiation field. These results indicate that the hot halo radiation fields potentially play an important role in quenching galaxies in cluster environments.

We study the stellar angular momentum of thousands of galaxies in the Illustris cosmological simulation, which captures gravitational and gas dynamics within galaxies, as well as feedback from stars and black holes. We find that the angular momentum of the simulated galaxies matches observations well, and in particular two distinct relations are found for late-type versus early-type galaxies. The relation for late-type galaxies corresponds to the value expected from full conservation of the specific angular momentum generated by cosmological tidal torques. The relation for early-type galaxies corresponds to retention of only ~30% of that, but we find that those early-type galaxies with low angular momentum at z=0 nevertheless reside at high redshift on the late-type relation. Some of them abruptly lose angular momentum during major mergers. To gain further insight, we explore the scaling relations in simulations where the galaxy formation physics is modified with respect to the fiducial model. We find that galactic winds with high mass-loading factors are essential for obtaining the high angular momentum relation typical for late-type galaxies, while AGN feedback largely operates in the opposite direction. Hence, feedback controls the stellar angular momentum of galaxies, and appears to be instrumental for establishing the Hubble sequence.

A method for improving the accuracy of hydrodynamical codes that use a moving Voronoi mesh is described. Our scheme is based on a new regularization scheme that constrains the mesh to be centroidal to high precision while still allowing the cells to move approximately with the local fluid velocity, thereby retaining the quasi-Lagrangian nature of the approach. Our regularization technique significantly reduces mesh noise that is attributed to changes in mesh topology and deviations from mesh regularity. We demonstrate the advantages of our method on various test problems, and note in particular improvements obtained in handling shear instabilities, mixing, and in angular momentum conservation. Calculations of adiabatic jets in which shear excites Kelvin Helmholtz instability show reduction of mesh noise and entropy generation. In contrast, simulations of the collapse and formation of an isolated disc galaxy are nearly unaffected, showing that numerical errors due to the choice of regularization do not impact the outcome in this case.

We have constructed merger trees for galaxies in the Illustris Simulation by directly tracking the baryonic content of subhalos. These merger trees are used to calculate the galaxy-galaxy merger rate as a function of descendant stellar mass, progenitor stellar mass ratio, and redshift. We demonstrate that the most appropriate definition for the mass ratio of a galaxy-galaxy merger consists in taking both progenitor masses at the time when the secondary progenitor reaches its maximum stellar mass. Additionally, we avoid effects from `orphaned' galaxies by allowing some objects to `skip' a snapshot when finding a descendant, and by only considering mergers which show a well-defined `infall' moment. Adopting these definitions, we obtain well-converged predictions for the galaxy-galaxy merger rate with the following main features, which are qualitatively similar to the halo-halo merger rate except for the last one: a strong correlation with redshift that evolves as $\sim (1+z)^{2.4-2.8}$, a power law with respect to mass ratio, and an increasing dependence on descendant stellar mass, which steepens significantly for descendant stellar masses greater than $\sim 2 \times 10^{11} \, {\rm M_{\odot}}$. These trends are consistent with observational constraints for medium-sized galaxies ($M_{\ast} \gtrsim 10^{10} \, {\rm M_{\odot}}$), but in tension with some recent observations of the close pair fraction for massive galaxies ($M_{\ast} \gtrsim 10^{11} \, {\rm M_{\odot}}$), which report a nearly constant or decreasing evolution with redshift. Finally, we provide a fitting function for the galaxy-galaxy merger rate which is accurate over a wide range of stellar masses, progenitor mass ratios, and redshifts.

We study the properties of gas in and around 10^12 solar mass halos at z=2 using a suite of high-resolution cosmological hydrodynamic 'zoom' simulations. We quantify the thermal and dynamical structure of these gaseous reservoirs in terms of their mean radial distributions and angular variability along different sightlines. With each halo simulated at three levels of increasing resolution, the highest reaching a baryon mass resolution of ~10,000 solar masses, we study the interaction of filamentary inflow and the quasi-static hot halo atmosphere. We highlight the discrepancy between the spatial resolution available in the halo gas as opposed to within the galaxy itself, and find that stream morphologies become increasingly complex at higher resolution, with large coherent flows revealing density and temperature structure at progressively smaller scales. Moreover, multiple gas components co-exist at the same radius within the halo, making radially averaged analyses misleading. This is particularly true where the hot, quasi-static, high entropy halo atmosphere interacts with cold, rapidly inflowing, low entropy accretion. We investigate the process of gas virialization and identify different regimes for the heating of gas as it accretes from the intergalactic medium. Haloes at this mass have a well-defined virial shock, associated with a sharp jump in temperature and entropy at ~1.25 r_vir. The presence, radius, and radial width of this boundary feature, however, vary not only from halo to halo, but also as a function of angular direction, covering roughly ~85% of the 4pi sphere. Our findings are relevant for the proper interpretation of observations pertaining to the circumgalactic medium, including evidence for large amounts of cold gas surrounding massive haloes at intermediate redshifts.

We investigate how the way galaxies acquire their gas across cosmic time in cosmological hydrodynamic simulations is modified by a comprehensive physical model for baryonic feedback processes. To do so, we compare two simulations -- with and without feedback -- both evolved with the moving mesh code AREPO. The feedback runs implement the full physics model of the Illustris simulation project, including star formation driven galactic winds and energetic feedback from supermassive blackholes. We explore: (a) the accretion rate of material contributing to the net growth of galaxies and originating directly from the intergalactic medium, finding that feedback strongly suppresses the raw, as well as the net, inflow of this "smooth mode" gas at all redshifts, regardless of the temperature history of newly acquired gas. (b) At the virial radius the temperature and radial flux of inflowing gas is largely unaffected at z=2. However, the spherical covering fraction of inflowing gas at 0.25 rvir decreases substantially, from more than 80% to less than 50%, while the rates of both inflow and outflow increase, indicative of recycling across this boundary. (c) The fractional contribution of smooth accretion to the total accretion rate is lower in the simulation with feedback, by roughly a factor of two across all redshifts. Moreover, the smooth component of gas with a cold temperature history, is entirely suppressed in the feedback run at z<1. (d) The amount of time taken by gas to cross from the virial radius to the galaxy -- the "halo transit time" -- increases in the presence of feedback by a factor of ~2-3, and is notably independent of halo mass. We discuss the possible implications of this invariance for theoretical models of hot halo gas cooling.

The presence of hot gaseous coronae around present-day massive spiral galaxies is a fundamental prediction of galaxy formation models. However, our observational knowledge remains scarce, since to date only four gaseous coronae were detected around spirals with massive stellar bodies ($\gtrsim2\times10^{11} \ \rm{M_{\odot}}$). To explore the hot coronae around lower mass spiral galaxies, we utilized Chandra X-ray observations of a sample of eight normal spiral galaxies with stellar masses of $(0.7-2.0)\times10^{11} \ \rm{M_{\odot}}$. Although statistically significant diffuse X-ray emission is not detected beyond the optical radii ($\sim20$ kpc) of the galaxies, we derive $3\sigma$ limits on the characteristics of the coronae. These limits, complemented with previous detections of NGC 1961 and NGC 6753, are used to probe the Illustris Simulation. The observed $3\sigma$ upper limits on the X-ray luminosities and gas masses exceed or are at the upper end of the model predictions. For NGC 1961 and NGC 6753 the observed gas temperatures, metal abundances, and electron density profiles broadly agree with those predicted by Illustris. These results hint that the physics modules of Illustris are broadly consistent with the observed properties of hot coronae around spiral galaxies. However, a shortcoming of Illustris is that massive black holes, mostly residing in giant ellipticals, give rise to powerful radio-mode AGN feedback, which results in under luminous coronae for ellipticals.

Massive, quiescent galaxies at high redshift have been found to be considerably more compact than galaxies of similar mass in the local universe. How these compact galaxies formed has yet to be determined, though several progenitor populations have been proposed. Here we investigate the formation processes and quantify the assembly histories of such galaxies in Illustris, a suite of hydrodynamical cosmological simulations encompassing a sufficiently large volume to include rare objects, while simultaneously resolving the internal structure of galaxies. We select massive (~10^11 solar masses) and compact (stellar half-mass radius < 2 kpc) galaxies from the simulation at z=2. Within the Illustris suite, we find that these quantities are not perfectly converged, but are reasonably reliable for our purposes. The resulting population is composed primarily of quiescent galaxies, but we also find several star-forming compact galaxies. The simulated compact galaxies are similar to observed galaxies in star formation activity and appearance. We follow their evolution at high redshift in the simulation and find that there are multiple pathways to form these compact galaxies, dominated by two mechanisms: (i) intense, centrally concentrated starbursts generally triggered by gas-rich major mergers between z~2-4, reducing the galaxies' half-mass radii by a factor of a few to below 2 kpc, and (ii) assembly at very early times when the universe was much denser; the galaxies formed compact and remained so until z~2.

Galactic feedback strongly affects the way galactic environments are enriched. We examine this connection by performing a suite of cosmological hydrodynamic simulations, exploring a range of parameters based on the galaxy formation model developed in Vogelsberger et al. 2013 (henceforth V13). We examine the effects of AGN feedback, wind mass loading, wind specific energy, and wind metal-loading on the properties of the circumgalactic medium (CGM) of galaxies with $M_\text{halo} > 10^{11} M_\odot$. Note that while the V13 model was tuned to match observations including the stellar mass function, no explicit tuning was done for the CGM. The wind energy per unit outflow mass has the most significant effect on the CGM enrichment. High energy winds launch metals far beyond the virial radius. AGN feedback also has a significant effect, but only at $z < 3$. We compare to high redshift HI and CIV observations. All our simulations produce the observed number of Damped Lyman-$\alpha$ Absorbers. At lower column density, several of our simulations produce enough Lyman Limit Systems (LLS) $100$ kpc from the galaxy, but in all cases the LLS abundance falls off with distance faster than observations, with too few LLS at $200$ kpc. Further, in all models the CIV abundance drops off too sharply with distance, with too little CIV $100$-$200$ kpc from the galaxy. Higher energy wind models produce more extended CIV but also produce less stars, in tension with star-formation rate density observations. This highlights the fact that circumgalactic observations are a strong constraint on galactic feedback models.

The formation and early evolution of circumstellar discs often occurs within dense, newborn stellar clusters. For the first time, we apply the moving-mesh code AREPO, to circumstellar discs in 3-D, focusing on disc-disc interactions that result from stellar fly-bys. Although a small fraction of stars are expected to undergo close approaches, the outcomes of the most violent encounters might leave an imprint on the discs and host stars that will influence both their orbits and their ability to form planets. We first construct well-behaved 3-D models of self-gravitating discs, and then create a suite of numerical experiments of parabolic encounters, exploring the effects of pericenter separation r_p, disc orientation and disc-star mass ratio (M_d/M_*) on the orbital evolution of the host stars. Close encounters (2r_p<~ disc radius) can truncate discs on very short time scales. If discs are massive, close encounters facilitate enough orbital angular momentum extraction to induce stellar capture. We find that for realistic primordial disc masses M_d<~0.1M_*, non-colliding encounters induce minor orbital changes, which is consistent with analytic calculations of encounters in the linear regime. The same disc masses produce entirely different results for grazing/colliding encounters. In the latter case, rapidly cooling discs lose orbital energy by radiating away the energy excess of the shock-heated gas, thus causing capture of the host stars into a bound orbit. In rare cases, a tight binary with a circumbinary disc forms as a result of this encounter.

We examine the kinematic structure of Damped Lyman-alpha Systems (DLAs) in a series of cosmological hydrodynamic simulations using the AREPO code. We are able to match the distribution of velocity widths of associated low ionisation metal absorbers substantially better than earlier work. Our simulations produce a population of DLAs dominated by halos with virial velocities around 70 km/s, consistent with a picture of relatively small, faint objects. In addition, we reproduce the observed correlation between velocity width and metallicity and the equivalent width distribution of SiII. Some discrepancies of moderate statistical significance remain; too many of our spectra show absorption concentrated at the edge of the profile and there are slight differences in the exact shape of the velocity width distribution. We show that the improvement over previous work is mostly due to our strong feedback from star formation and our detailed modelling of the metal ionisation state.

Understanding the physical processes that drive star formation is a key challenge for galaxy formation models. In this article we study the tight correlation between the star formation rate (SFR) and stellar mass of galaxies at a given redshift, how halo growth influences star formation, and star formation histories of individual galaxies. We study these topics using Illustris, a state-of-the-art cosmological hydrodynamical simulation of galaxy formation. Illustris reproduces the observed relation (the star formation main sequence; SFMS) between SFR and stellar mass at redshifts z=0 and z=4, but at intermediate redshifts of z~2, the simulated SFMS has a significantly lower normalisation than reported by observations. The scatter in the relation is consistent with the observed scatter. However, the fraction of outliers above the SFR-stellar mass relation in Illustris is less than that observed. Galaxies with halo masses of ~10^{12} solar masses dominate the SFR density of the Universe, in agreement with the results of abundance matching. Furthermore, more-massive galaxies tend to form the bulk of their stars at high redshift, which indicates that `downsizing' occurs in Illustris. We also studied the star formation histories of individual galaxies, including the use of a principal component analysis decomposition. We find that for fixed stellar mass, galaxies that form earlier have more-massive black holes at z=0, indicating that star formation and black hole growth are tightly linked processes in Illustris. While many of the properties of normal star-forming galaxies are well-reproduced in the Illustris simulation, forming a realistic population of starbursts will likely require higher resolution and probably a more sophisticated treatment of star formation and feedback from stars and black holes.

We present our methods for generating a catalog of 7,000 synthetic images and 40,000 integrated spectra of redshift z = 0 galaxies from the Illustris Simulation. The mock data products are produced by using stellar population synthesis models to assign spectral energy distributions (SED) to each star particle in the galaxies. The resulting synthetic images and integrated SEDs therefore properly reflect the spatial distribution, stellar metallicity distribution, and star formation history of the galaxies. From the synthetic data products it is possible to produce monochromatic or color-composite images, perform SED fitting, classify morphology, determine galaxy structural properties, and evaluate the impacts of galaxy viewing angle. The main contribution of this paper is to describe the production, format, and composition of the image catalog that makes up the Illustris Simulation Obsevatory. As a demonstration of this resource, we derive galactic stellar mass estimates by applying the SED fitting code FAST to the synthetic galaxy products, and compare the derived stellar masses against the true stellar masses from the simulation. We find from this idealized experiment that systematic biases exist in the photometrically derived stellar mass values that can be reduced by using a fixed metallicity in conjunction with a minimum galaxy age restriction.

Observationally, the fraction of blue satellite galaxies decreases steeply with host halo mass, and their radial distribution around central galaxies is significantly shallower in massive (M_* >10e11M_sun) than in Milky Way like systems. Theoretical models, based primarily on semi-analytical techniques, have had a long-standing problem with reproducing these trends, instead predicting too few blue satellites in general but also estimating a radial distribution that is too shallow, regardless of primary mass. In this Letter, we use the Illustris cosmological simulation to study the properties of satellite galaxies around isolated primaries. For the first time, we find good agreement between theory and observations. We identify the main source of this success relative to earlier work to be a consequence of the large gas contents of satellites at infall, a factor ~5-10 times larger than in semi-analytical models. Because of their relatively large gas reservoirs, satellites can continue to form stars long after infall, with a typical timescale for star-formation to be quenched ~2 Gyr in groups but more than ~5 Gyr for satellites around Milky Way like primaries. The gas contents we infer are consistent with z=0 observations of HI gas in galaxies, although we find large discrepancies among reported values in the literature. A testable prediction of our model is that the gas-to-stellar mass ratio of satellite progenitors should vary only weakly with cosmic time.

We introduce the Illustris Project, a series of large-scale hydrodynamical simulations of galaxy formation. The highest resolution simulation, Illustris-1, covers a volume of $(106.5\,{\rm Mpc})^3$, has a dark mass resolution of ${6.26 \times 10^{6}\,{\rm M}_\odot}$, and an initial baryonic matter mass resolution of ${1.26 \times 10^{6}\,{\rm M}_\odot}$. At $z=0$ gravitational forces are softened on scales of $710\,{\rm pc}$, and the smallest hydrodynamical gas cells have an extent of $48\,{\rm pc}$. We follow the dynamical evolution of $2\times 1820^3$ resolution elements and in addition passively evolve $1820^3$ Monte Carlo tracer particles reaching a total particle count of more than $18$ billion. The galaxy formation model includes: primordial and metal-line cooling with self-shielding corrections, stellar evolution, stellar feedback, gas recycling, chemical enrichment, supermassive black hole growth, and feedback from active galactic nuclei. At $z=0$ our simulation volume contains about $40,000$ well-resolved galaxies covering a diverse range of morphologies and colours including early-type, late-type and irregular galaxies. The simulation reproduces reasonably well the cosmic star formation rate density, the galaxy luminosity function, and baryon conversion efficiency at $z=0$. It also qualitatively captures the impact of galaxy environment on the red fractions of galaxies. The internal velocity structure of selected well-resolved disk galaxies obeys the stellar and baryonic Tully-Fisher relation together with flat circular velocity curves. In the well-resolved regime the simulation reproduces the observed mix of early-type and late-type galaxies. Our model predicts a halo mass dependent impact of baryonic effects on the halo mass function and the masses of haloes caused by feedback from supernova and active galactic nuclei.

We present the first cosmological simulations of dwarf galaxies, which include dark matter self-interactions and baryons. We study two dwarf galaxies within cold dark matter, and four different elastic self-interacting scenarios with constant and velocity-dependent cross sections, motivated by a new force in the hidden dark matter sector. Our highest resolution simulation has a baryonic mass resolution of $1.8\times 10^2\,{\rm M}_\odot$ and a gravitational softening length of $34\,{\rm pc}$ at $z=0$. In this first study we focus on the regime of mostly isolated dwarf galaxies with halo masses $\sim10^{10}\,{\rm M}_\odot$ where dark matter dynamically dominates even at sub-kpc scales. We find that while the global properties of galaxies of this scale are minimally affected by allowed self-interactions, their internal structures change significantly if the cross section is large enough within the inner sub-kpc region. In these dark-matter-dominated systems, self-scattering ties the shape of the stellar distribution to that of the dark matter distribution. In particular, we find that the stellar core radius is closely related to the dark matter core radius generated by self-interactions. Dark matter collisions lead to dwarf galaxies with larger stellar cores and smaller stellar central densities compared to the cold dark matter case. The central metallicity within $1\,{\rm kpc}$ is also larger by up to $\sim 15\%$ in the former case. We conclude that the mass distribution, and characteristics of the central stars in dwarf galaxies can potentially be used to probe the self-interacting nature of dark matter.

We examine the abundance, clustering and metallicity of Damped Lyman-alpha Absorbers (DLAs) in a suite of hydrodynamic cosmological simulations using the moving mesh code AREPO. We incorporate models of supernova and AGN feedback, as well as molecular hydrogen formation. We compare our simulations to the column density distribution function at $z=3$, the total DLA abundance at $z=2-4$, the measured DLA bias at $z=2.3$ and the DLA metallicity distribution at $z=2-4$. Our preferred models produce populations of DLAs in good agreement with most of these observations. The exception is the DLA abundance at $z < 3$, which we show requires stronger feedback in $10^{11-12} \, h^{-1} M_\odot$ mass halos. While the DLA population probes a wide range of halo masses, we find the cross-section is dominated by halos of mass $10^{10} - 10^{11} \, h^{-1} M_\odot$ and virial velocities $50 - 100 \;\mathrm{km/s}$. The simulated DLA population has a linear theory bias of $1.7$, whereas the observations require $2.17 \pm 0.2$. We show that non-linear growth increases the bias in our simulations to $2.3$ at $k=1\; \mathrm{Mpc/}h$, the smallest scale observed. The scale-dependence of the bias is, however, very different in the simulations compared against the observations. We show that, of the observations we consider, the DLA abundance and column density function provide the strongest constraints on the feedback model.

We use the Illustris Simulations to gain insight into the build-up of the outer, low-surface brightness regions which surround galaxies. We characterize the stellar haloes by means of the logarithmic slope of the spherically-averaged stellar density profiles, alphaSTARS at z=0, and we relate these slopes to the properties of the underlying Dark-Matter (DM) haloes, their central galaxies, and their assembly histories. We analyze a sample of ~5,000 galaxies resolved with more than 5x10^4 particles each, and spanning a variety of morphologies and halo masses (3x10^11 < Mvir < 10^14 Msun). We find a strong trend between stellar halo slope and total halo mass, where more massive objects have shallower stellar haloes than the less massive ones (-5.5 \pm 0.5 < alphaSTARS <-3.5 \pm 0.2 in the studied mass range). At fixed halo mass, we show that disk-like, blue, young, and more massive galaxies are surrounded by significantly steeper stellar haloes than elliptical, red, older, and less massive galaxies. Overall, the stellar density profiles fall off much more steeply than the underlying DM, and no clear trend holds between stellar slope and DM halo concentration. However, DM haloes which formed more recently, or which accreted larger fractions of stellar mass from infalling satellites, exhibit shallower stellar haloes than their older analogs with similar masses, by up to Delta(alphaSTARS) ~ 0.5-0.7. Our findings, combined with the most recent measurements of the strikingly different stellar power-law indexes for M31 and the Milky Way, appear to favour a massive M31, and a Milky Way characterized by a much quieter accretion history over the past 10 Gyrs than its companion.

We present an overview of galaxy evolution across cosmic time in the Illustris Simulation. Illustris is an N-body/hydrodynamical simulation that evolves 2*1820^3 resolution elements in a (106.5Mpc)^3 box from cosmological initial conditions down to z=0 using the AREPO moving-mesh code. The simulation uses a state-of-the-art set of physical models for galaxy formation that was tuned to reproduce the z=0 stellar mass function and the history of the cosmic star-formation rate density. We find that Illustris successfully reproduces a plethora of observations of galaxy populations at various redshifts, for which no tuning was performed, and provide predictions for future observations. In particular, we discuss (a) the buildup of galactic mass, showing stellar mass functions and the relations between stellar mass and halo mass from z=7 to z=0, (b) galaxy number density profiles around massive central galaxies out to z=4, (c) the gas and total baryon content of both galaxies and their halos for different redshifts, and as a function of mass and radius, and (d) the evolution of galaxy specific star-formation rates up to z=8. In addition, we (i) present a qualitative analysis of galaxy morphologies from z=5 to z=0, for the stellar as well as the gaseous components, and their appearance in HST mock observations, (ii) follow galaxies selected at z=2 to their z=0 descendants, and quantify their growth and merger histories, and (iii) track massive z=0 galaxies to high redshift and study their joint evolution in star-formation activity and compactness. We conclude with a discussion of several disagreements with observations, and lay out possible directions for future research.

We consider three extensions of the Navarro, Frenk and White (NFW) profile and investigate the intrinsic degeneracies among the density profile parameters on the gravitational lensing effect of satellite galaxies on highly magnified Einstein rings. In particular, we find that the gravitational imaging technique can be used to exclude specific regions of the considered parameter space, and therefore, models that predict a large number of satellites in those regions. By comparing the lensing degeneracy with the intrinsic density profile degeneracies, we show that theoretical predictions based on fits that are dominated by the density profile at larger radii may significantly over- or underestimate the number of satellites that are detectable with gravitational lensing. Finally, using the previously reported detection of a satellite in the gravitational lens system JVAS B1938+666 as an example, we derive for this detected satellite values of r_max and v_max that are, for each considered profile, consistent within 1sigma with the parameters found for the luminous dwarf satellites of the Milky Way and with a mass density slope gamma < 1.6. We also find that the mass of the satellite within the Einstein radius as measured using gravitational lensing is stable against assumptions on the substructure profile. In the future thanks to the increased angular resolution of very long baseline interferometry at radio wavelengths and of the E-ELT in the optical we will be able to set tighter constraints on the number of allowed substructure profiles.

Galaxy mergers have been investigated for decades using smoothed particle hydrodynamics (SPH), but recent work highlighting inaccuracies inherent in the traditional SPH technique calls into question the reliability of previous studies. We explore this issue by comparing a suite of Gadget-3 SPH simulations of idealised (i.e., non-cosmological) isolated discs and galaxy mergers with otherwise identical calculations performed using the moving-mesh code Arepo. When black hole (BH) accretion and active galactic nucleus (AGN) feedback are not included, the star formation histories (SFHs) obtained from the two codes agree well. When BHs are included, the code- and resolution-dependent variations in the SFHs are more significant, but the agreement is still good, and the stellar mass formed over the course of a simulation is robust to variations in the numerical method. During a merger, the gas morphology and phase structure are initially similar prior to the starburst phase. However, once a hot gaseous halo has formed from shock heating and AGN feedback (when included), the agreement is less good. In particular, during the post-starburst phase, the SPH simulations feature more prominent hot gaseous haloes and spurious clumps, whereas with Arepo, gas clumps and filaments are less apparent and the hot halo gas can cool more efficiently. We discuss the origin of these differences and explain why the SPH technique yields trustworthy results for some applications (such as the idealised isolated disc and galaxy merger simulations presented here) but not others (e.g., gas flows onto galaxies in cosmological hydrodynamical simulations).

Magnetic fields play an important role in many astrophysical systems and a detailed understanding of their impact on the gas dynamics requires robust numerical simulations. Here we present a new method to evolve the ideal magnetohydrodynamic (MHD) equations on unstructured static and moving meshes that preserves the magnetic field divergence-free constraint to machine precision. The method overcomes the major problems of using a cleaning scheme on the magnetic fields instead, which is non-conservative, not fully Galilean invariant, does not eliminate divergence errors completely, and may produce incorrect jumps across shocks. Our new method is a generalization of the constrained transport (CT) algorithm used to enforce the $\nabla\cdot \mathbf{B}=0$ condition on fixed Cartesian grids. Preserving $\nabla\cdot \mathbf{B}=0$ at the discretized level is necessary to maintain the orthogonality between the Lorentz force and $\mathbf{B}$. The possibility of performing CT on a moving mesh provides several advantages over static mesh methods due to the quasi-Lagrangian nature of the former (i.e., the mesh generating points move with the flow), such as making the simulation automatically adaptive and significantly reducing advection errors. Our method preserves magnetic fields and fluid quantities in pure advection exactly.

We demonstrate that self-interacting dark matter models with interactions mediated by light particles can have significant deviations in the matter power-spectrum and detailed structure of galactic halos when compared to a standard cold dark matter scenario. While these deviations can take the form of suppression of small scale structure that are in some ways similar to that of warm dark matter, the self-interacting models have a much wider range of possible phenomenology. A long-range force in the dark matter can introduce multiple scales to the initial power spectrum, in the form of dark acoustic oscillations and an exponential cut-off in the power spectrum. Using simulations we show that the impact of these scales can remain observationally relevant up to the present day. Furthermore, the self-interaction can continue to modify the small-scale structure of the dark matter halos, reducing their central densities and creating a dark matter core. The resulting phenomenology is unique to this type of models.

Previous simulations of the growth of cosmic structures have broadly reproduced the 'cosmic web' of galaxies that we see in the Universe, but failed to create a mixed population of elliptical and spiral galaxies due to numerical inaccuracies and incomplete physical models. Moreover, because of computational constraints, they were unable to track the small scale evolution of gas and stars to the present epoch within a representative portion of the Universe. Here we report a simulation that starts 12 million years after the Big Bang, and traces 13 billion years of cosmic evolution with 12 billion resolution elements in a volume of $(106.5\,{\rm Mpc})^3$. It yields a reasonable population of ellipticals and spirals, reproduces the distribution of galaxies in clusters and statistics of hydrogen on large scales, and at the same time the metal and hydrogen content of galaxies on small scales.

We investigate how different cosmological parameters, such as those delivered by the WMAP and Planck missions, affect the nature and evolution of dark matter halo substructure. We use a series of flat $\Lambda$ cold dark matter ($\Lambda$CDM) cosmological $N$-body simulations of structure formation, each with a different power spectrum but the same initial white noise field. Our fiducial simulation is based on parameters from the WMAP 7th year cosmology. We then systematically vary the spectral index, $n_s$, matter density, $\Omega_M$, and normalization of the power spectrum, $\sigma_8$, for 7 unique simulations. Across these, we study variations in the subhalo mass function, mass fraction, maximum circular velocity function, spatial distribution, concentration, formation times, accretion times, and peak mass. We eliminate dependence of subhalo properties on host halo mass and average over many hosts to reduce variance. While the "same" subhalos from identical initial overdensity peaks in higher $\sigma_8, n_s$, and $\Omega_m$ simulations accrete earlier and end up less massive and closer to the halo center at $z=0$, the process of continuous subhalo accretion and destruction leads to a steady state distribution of these properties across all subhalos in a given host. This steady state mechanism eliminates cosmological dependence on all properties listed above except subhalo concentration and $V_{max}$, which remain greater for higher $\sigma_8, n_s$ and $\Omega_m$ simulations, and subhalo formation time, which remains earlier. We also find that the numerical technique for computing scale radius and the halo finder used can significantly affect the concentration-mass relationship computed for a simulation.

Dynamic Complexity is a phenomenon exhibited by a nonlinearly interacting system within which multitudes of different sizes of large scale coherent structures emerge, resulting in a globally nonlinear stochastic behavior vastly different from that could be surmised from the underlying equations of interaction. The hallmark of such nonlinear, complex phenomena is the appearance of intermittent fluctuating events with the mixing and distributions of correlated structures at all scales. We briefly review here a relatively recent method, ROMA (rank-ordered multifractal analysis), explicitly constructed to analyze the intricate details of the distribution and scaling of such types of intermittent structures. This method is then applied to the analyses of selected examples related to the dynamical plasmas of the cusp region of the magnetosphere, velocity fluctuations of classical hydrodynamic turbulence, and the distribution of the structures of the cosmic gas obtained through large scale, moving mesh simulations. Differences and similarities of the analyzed results among these complex systems will be contrasted and highlighted. The first two examples have direct relevance to the geospace environment and are summaries of previously reported findings. The third example on the cosmic gas, though involving phenomena much larger in spatiotemporal scales, with its highly compressible turbulent behavior and the unique simulation technique employed in generating the data, provides direct motivations of applying such analysis to studies of similar multifractal processes in various extreme environments. These new results are both exciting and intriguing.

We present a new comprehensive model of the physics of galaxy formation designed for large-scale hydrodynamical simulations of structure formation using the moving mesh code AREPO. Our model includes primordial and metal line cooling with self-shielding corrections, stellar evolution and feedback processes, gas recycling, chemical enrichment, a novel subgrid model for the metal loading of outflows, black hole (BH) seeding, BH growth and merging procedures, quasar- and radio-mode feedback, and a prescription for radiative electro-magnetic (EM) feedback from active galactic nuclei (AGN). The metal mass loading of outflows can be adjusted independently of the wind mass loading. This is required to simultaneously reproduce the stellar mass content of low mass haloes and their gas oxygen abundances. Radiative EM AGN feedback is implemented assuming an average spectral energy distribution and a luminosity-dependent scaling of obscuration effects. This form of feedback suppresses star formation more efficiently than continuous thermal quasar-mode feedback alone, but is less efficient than mechanical radio-mode feedback in regulating star formation in massive haloes. We contrast simulation predictions for different variants of our galaxy formation model with key observations. Our best match model reproduces, among other things, the cosmic star formation history, the stellar mass function, the stellar mass - halo mass relation, g-, r-, i-, z-band SDSS galaxy luminosity functions, and the Tully-Fisher relation. We can achieve this success only if we invoke very strong forms of stellar and AGN feedback such that star formation is adequately reduced in both low and high mass systems. In particular, the strength of radio-mode feedback needs to be increased significantly compared to previous studies to suppress efficient cooling in massive, metal-enriched haloes.

Galactic winds are a key physical mechanism for understanding galaxy formation and evolution, yet empirical and theoretical constraints for the character of winds are limited and discrepant. Recent empirical models find that local star-forming galaxies have a deficit of oxygen that scales with galaxy stellar mass. The oxygen deficit provides unique empirical constraints on the magnitude of mass loss, composition of outflowing material and metal reaccretion onto galaxies. We formulate the oxygen deficit constraints so they may be easily implemented into theoretical models of galaxy evolution. We parameterize an effective metal loading factor which combines the uncertainties of metal outflows and metal reaccretion into a single function of galaxy virial velocity. We determine the effective metal loading factor by forward-fitting the oxygen deficit. The effective metal loading factor we derive has important implications for the implementation of mass loss in models of galaxy evolution.

We present two numerical schemes for passive tracer particles in the hydrodynamical moving-mesh code AREPO, and compare their performance for various problems, from simple setups to cosmological simulations. The purpose of tracer particles is to allow the flow to be followed in a Lagrangian way, tracing the evolution of the fluid with time, and allowing the thermodynamical history of individual fluid parcels to be recorded. We find that the commonly-used `velocity field tracers', which are advected using the fluid velocity field, do not in general follow the mass flow correctly, and explain why this is the case. This method can result in orders-of-magnitude biases in simulations of driven turbulence and in cosmological simulations, rendering the velocity field tracers inappropriate for following these flows. We then discuss a novel implementation of `Monte Carlo tracers', which are moved along with fluid cells, and are exchanged probabilistically between them following the mass flux. This method reproduces the mass distribution of the fluid correctly. The main limitation of this approach is that it is more diffusive than the fluid itself. Nonetheless, we show that this novel approach is more reliable than what has been employed previously and demonstrate that it is appropriate for following hydrodynamical flows in mesh-based codes. The Monte Carlo tracers can also naturally be transferred between fluid cells and other types of particles, such as stellar particles, so that the mass flow in cosmological simulations can be followed in its entirety.

Luminous X-ray gas coronae in the dark matter halos of massive spiral galaxies are a fundamental prediction of structure formation models, yet only a few such coronae have been detected so far. In this paper, we study the hot X-ray coronae beyond the optical disks of two normal massive spirals, NGC1961 and NGC6753. Based on XMM-Newton X-ray observations, hot gaseous emission is detected to ~60 kpc - well beyond their optical radii. The hot gas has a best-fit temperature of kT~0.6 keV and an abundance of ~0.1 Solar, and exhibits a fairly uniform distribution, suggesting that the quasi-static gas resides in hydrostatic equilibrium in the potential well of the galaxies. The bolometric luminosity of the gas in the (0.05-0.15)r_200 region (r_200 is the virial radius) is ~6e40 erg/s for both galaxies. The baryon mass fractions of NGC1961 and NGC6753 are f_b~0.1, which fall short of the cosmic baryon fraction. The hot coronae around NGC1961 and NGC6753 offer an excellent basis to probe structure formation simulations. To this end, the observations are confronted with the moving mesh code Arepo and the smoothed particle hydrodynamics code Gadget. Although neither model gives a perfect description, the observed luminosities, gas masses, and abundances favor the Arepo code. Moreover, the shape and the normalization of the observed density profiles are better reproduced by Arepo within ~0.5r_200. However, neither model incorporates efficient feedback from supermassive black holes or supernovae, which could alter the simulated properties of the X-ray coronae. With the further advance of numerical models, the present observations will be essential in constraining the feedback effects in structure formation simulations.

A discontinuous Galerkin (DG) method suitable for large-scale astrophysical simulations on Cartesian meshes as well as arbitrary static and moving Voronoi meshes is presented. Most major astrophysical fluid dynamics codes use a finite volume (FV) approach. We demonstrate that the DG technique offers distinct advantages over FV formulations on both static and moving meshes. The DG method is also easily generalized to higher than second-order accuracy without requiring the use of extended stencils to estimate derivatives (thereby making the scheme highly parallelizable). We implement the technique in the AREPO code for solving the fluid and the magnetohydrodynamic (MHD) equations. By examining various test problems, we show that our new formulation provides improved accuracy over FV approaches of the same order, and reduces post-shock oscillations and artificial diffusion of angular momentum. In addition, the DG method makes it possible to represent magnetic fields in a locally divergence-free way, improving the stability of MHD simulations and moderating global divergence errors, and is a viable alternative for solving the MHD equations on meshes where Constrained-Transport (CT) cannot be applied. We find that the DG procedure on a moving mesh is more sensitive to the choice of slope limiter than is its FV method counterpart. Therefore, future work to improve the performance of the DG scheme even further will likely involve the design of optimal slope limiters. As presently constructed, our technique offers the potential of improved accuracy in astrophysical simulations using the moving mesh AREPO code as well as those employing adaptive mesh refinement (AMR).

We present a multi-epoch analysis of the galaxy populations formed within the cosmological hydrodynamical simulations presented in Vogelsberger et al. (2013). These simulations explore the performance of a recently implemented feedback model which includes primordial and metal line radiative cooling with self-shielding corrections; stellar evolution with associated mass loss and chemical enrichment; feedback by stellar winds; black hole seeding, growth and merging; and AGN quasar- and radio-mode heating with a phenomenological prescription for AGN electro-magnetic feedback. We illustrate the impact of the model parameter choices on the resulting simulated galaxy population properties at high and intermediate redshifts. We demonstrate that our scheme is capable of producing galaxy populations that broadly reproduce the observed galaxy stellar mass function extending from redshift z=0 to z=3. We also characterise the evolving galactic B-band luminosity function, stellar mass to halo mass ratio, star formation main sequence, Tully-Fisher relation, and gas-phase mass-metallicity relation and confront them against recent observational estimates. This detailed comparison allows us to validate elements of our feedback model, while also identifying areas of tension that will be addressed in future work.

We quantify the impact of massive neutrinos on the statistics of low density regions in the intergalactic medium (IGM) as probed by the Lyman-alpha forest at redshifts z=2.2--4. Based on mock but realistic quasar (QSO) spectra extracted from hydrodynamic simulations with cold dark matter, baryons and neutrinos, we find that the probability distribution of weak Lyman-alpha absorption features, as sampled by Lyman-alpha flux regions at high transmissivity, is strongly affected by the presence of massive neutrinos. We show that systematic errors affecting the Lyman-alpha forest reduce but do not erase the neutrino signal. Using the Fisher matrix formalism, we conclude that the sum of the neutrino masses can be measured, using the method proposed in this paper, with a precision smaller than 0.4 eV using a catalog of 200 high resolution (S/N~100) QSO spectra. This number reduces to 0.27 eV by making use of reasonable priors in the other parameters that also affect the statistics of the high transitivity regions of the Lyman-alpha forest. The constraints obtained with this method can be combined with independent bounds from the CMB, large scale structures and measurements of the matter power spectrum from the Lyman-alpha forest to produce tighter upper limits on the sum of the masses of the neutrinos.

We investigate the nature of gas accretion onto haloes and galaxies at z=2 using cosmological hydrodynamic simulations run with the moving mesh code AREPO. Implementing a Monte Carlo tracer particle scheme to determine the origin and thermodynamic history of accreting gas, we make quantitative comparisons to an otherwise identical simulation run with the smoothed particle hydrodynamics (SPH) code GADGET-3. Contrasting these two numerical approaches, we find significant physical differences in the thermodynamic history of accreted gas in haloes above 10^10.5 solar masses. In agreement with previous work, GADGET simulations show a cold fraction near unity for galaxies forming in massive haloes, implying that only a small percentage of accreted gas heats to an appreciable fraction of the virial temperature during accretion. The same galaxies in AREPO show a much lower cold fraction, <20% in haloes above 10^11 solar masses. This results from a hot gas accretion rate which, at this same halo mass, is an order of magnitude larger than with GADGET, while the cold accretion rate is also lower. These discrepancies increase for more massive systems, and we explain both as due to numerical inaccuracies in the standard formulation of SPH. We also observe that the relatively sharp transition from cold to hot mode dominated accretion, at a halo mass of ~10^11, is a consequence of comparing past gas temperatures to a constant threshold value independent of virial temperature. Examining the spatial distribution of accreting gas, we find that gas filaments in GADGET tend to remain collimated and flow coherently to small radii, or artificially fragment and form a large number of purely numerical "blobs". Similar gas streams in AREPO show increased heating and disruption at 0.25-0.5 virial radii and contribute to the hot gas accretion rate in a manner distinct from classical cooling flows.

The causes of spiral structure in galaxies remain uncertain. Leaving aside the grand bisymmetric spirals with their own well-known complications, here we consider the possibility that multi-armed spiral features originate from density inhomogeneities orbiting within disks. Using high-resolution N-body simulations, we follow the motions of stars under the influence of gravity, and show that mass concentrations with properties similar to those of giant molecular clouds can induce the development of spiral arms through a process termed swing amplification. However, unlike in earlier work, we demonstrate that the eventual response of the disk can be highly non-linear, significantly modifying the formation and longevity of the resulting patterns. Contrary to expectations, ragged spiral structures can thus survive at least in a statistical sense long after the original perturbing influence has been removed.

Self-Interacting Dark Matter is an attractive alternative to the Cold Dark Matter paradigm only if it is able to substantially reduce the central densities of dwarf-size haloes while keeping the densities and shapes of cluster-size haloes within current constraints. Given the seemingly stringent nature of the latter, it was thought for nearly a decade that SIDM would be viable only if the cross section for self-scattering was strongly velocity-dependent. However, it has recently been suggested that a constant cross section per unit mass of sigma_T/m~0.1cm^2/g is sufficient to accomplish the desired effect. We explicitly investigate this claim using high resolution cosmological simulations of a Milky-Way size halo and find that, similarly to the Cold Dark Matter case, such cross section produces a population of massive subhaloes that is inconsistent with the kinematics of the classical dwarf spheroidals, in particular with the inferred slopes of the mass profiles of Fornax and Sculptor. This problem is resolved if sigma_T/m~1cm^2/g at the dwarf spheroidal scales. Since this value is likely inconsistent with the halo shapes of several clusters, our results leave only a small window open for a velocity-independent Self-Interacting Dark Matter model to work as a distinct alternative to Cold Dark Matter.

Self-interacting dark matter offers an interesting alternative to collisionless dark matter because of its ability to preserve the large-scale success of the cold dark matter model, while seemingly solving its challenges on small scales. We present here the first study of the expected dark matter detection signal taking into account different self-scattering models. We demonstrate that models with constant and velocity dependent cross sections, which are consistent with observational constraints, lead to distinct signatures in the velocity distribution, because non-thermalised features found in the cold dark matter distribution are thermalised through particle scattering. Depending on the model, self-interaction can lead to a 10% reduction of the recoil rates at high energies, corresponding to a minimum speed that can cause recoil larger than 300 km/s, compared to the cold dark matter case. At lower energies these differences are smaller than 5% for all models. The amplitude of the annual modulation signal can increase by up to 25%, and the day of maximum amplitude can shift by about two weeks with respect to the cold dark matter expectation. Furthermore, the exact day of phase reversal of the modulation signal can also differ by about a week between the different models. In general, models with velocity dependent cross sections peaking at the typical velocities of dwarf galaxies lead only to minor changes in the detection signals, whereas allowed constant cross section models lead to significant changes. We conclude that different self-interacting dark matter scenarios might be distinguished from each other through the details of direct detection signals. Furthermore, detailed constraints on the intrinsic properties of dark matter based on null detections, should take into account the possibility of self-scattering and the resulting effects on the detector signal.

We examine the distribution of neutral hydrogen in cosmological simulations carried out with the new moving-mesh code AREPO and compare it with the corresponding GADGET simulations based on the smoothed particle hydrodynamics (SPH) technique. The two codes use identical gravity solvers and baryonic physics implementations, but very different methods for solving the Euler equations, allowing us to assess how numerical effects associated with the hydro-solver impact the results of simulations. Here we focus on an analysis of the neutral gas, as detected in quasar absorption lines. We find that the high column density regime probed by Damped Lyman-alpha (DLA) and Lyman Limit Systems (LLS) exhibits significant differences between the codes. GADGET produces spurious artefacts in large halos in the form of gaseous clumps, boosting the LLS cross-section. Furthermore, it forms halos with denser central baryonic cores than AREPO, which leads to a substantially greater DLA cross-section from smaller halos. AREPO thus produces a significantly lower cumulative abundance of DLAs, which is intriguingly in much closer agreement with observations. The column density function, however, is not altered enough to significantly reduce the discrepancy with the observed value. For the low column density gas probed by the Lyman-alpha forest, the codes differ only at the level of a few percent, suggesting that this regime is quite well described by both methods, a fact that is reassuring for the many Lyman-alpha studies carried out with SPH thus far. While the residual differences are smaller than the errors on current Lyman-alpha forest data, we note that this will likely change for future precision experiments.

We present a review of the current state of the art of cosmological dark matter simulations, with particular emphasis on the implications for dark matter detection efforts and studies of dark energy. This review is intended both for particle physicists, who may find the cosmological simulation literature opaque or confusing, and for astro-physicists, who may not be familiar with the role of simulations for observational and experimental probes of dark matter and dark energy. Our work is complementary to the contribution by M. Baldi in this issue, which focuses on the treatment of dark energy and cosmic acceleration in dedicated N-body simulations. Truly massive dark matter-only simulations are being conducted on national supercomputing centers, employing from several billion to over half a trillion particles to simulate the formation and evolution of cosmologically representative volumes (cosmic scale) or to zoom in on individual halos (cluster and galactic scale). These simulations cost millions of core-hours, require tens to hundreds of terabytes of memory, and use up to petabytes of disk storage. The field is quite internationally diverse, with top simulations having been run in China, France, Germany, Korea, Spain, and the USA. Predictions from such simulations touch on almost every aspect of dark matter and dark energy studies, and we give a comprehensive overview of this connection. We also discuss the limitations of the cold and collisionless DM-only approach, and describe in some detail efforts to include different particle physics as well as baryonic physics in cosmological galaxy formation simulations, including a discussion of recent results highlighting how the distribution of dark matter in halos may be altered. We end with an outlook for the next decade, presenting our view of how the field can be expected to progress. (abridged)

We compare the structural properties of galaxies formed in cosmological simulations using the smoothed particle hydrodynamics (SPH) code GADGET with those using the moving-mesh code AREPO. Both codes employ identical gravity solvers and the same sub-resolution physics but use very different methods to track the hydrodynamic evolution of gas. This permits us to isolate the effects of the hydro solver on the formation and evolution of galactic gas disks in GADGET and AREPO haloes with comparable numerical resolution. In a matching sample of GADGET and AREPO haloes we fit simulated gas disks with exponential profiles. We find that the cold gas disks formed using the moving mesh approach have systematically larger disk scale lengths and higher specific angular momenta than their GADGET counterparts across a wide range in halo masses. For low mass galaxies differences between the properties of the simulated galaxy disks are caused by an insufficient number of resolution elements which lead to the artificial angular momentum transfer in our SPH calculation. We however find that galactic disks formed in massive halos, resolved with 10^6 particles/cells, are still systematically smaller in the GADGET run by a factor of ~2. The reasons for this are: 1) The excessive heating of haloes close to the cooling radius due to spurious dissipation of the subsonic turbulence in GADGET; and 2) The efficient delivery of low angular momentum gaseous blobs to the bottom of the potential well. While this large population of gaseous blobs in GADGET originates from the filaments which are pressure confined and fragment due to the SPH surface tension while infalling into hot halo atmospheres, it is essentially absent in the moving mesh calculation, clearly indicating numerical rather than physical origin of the blob material.

We discuss cosmological hydrodynamic simulations of galaxy formation performed with the new moving-mesh code AREPO, which promises higher accuracy compared with the traditional SPH technique that has been widely employed for this problem. We use an identical set of physics in corresponding simulations carried out with the well-tested SPH code GADGET, adopting also the same high-resolution gravity solver. We are thus able to compare both simulation sets on an object-by-object basis, allowing us to cleanly isolate the impact of different hydrodynamical methods on galaxy and halo properties. In accompanying papers, we focus on an analysis of the global baryonic statistics predicted by the simulation codes, (Vogelsberger et al. 2011) and complementary idealized simulations that highlight the differences between the hydrodynamical schemes (Sijacki et al. 2011). Here we investigate their influence on the baryonic properties of simulated galaxies and their surrounding haloes. We find that AREPO leads to significantly higher star formation rates for galaxies in massive haloes and to more extended gaseous disks in galaxies, which also feature a thinner and smoother morphology than their GADGET counterparts. Consequently, galaxies formed in AREPO have larger sizes and higher specific angular momentum than their SPH correspondents. The more efficient cooling flows in AREPO yield higher densities and lower entropies in halo centers (and the opposite trend in halo outskirts) leading to higher star formation rates of massive galaxies. While both codes agree to acceptable accuracy on a number of baryonic properties of cosmic structures, our results clearly demonstrate that galaxy formation simulations greatly benefit from the use of more accurate hydrodynamical techniques such as AREPO.

The Fermi-LAT collaboration has recently reported the detection of angular power above the photon noise level in the diffuse gamma-ray background between 1 and 50 GeV. Such signal can be used to constrain a possible contribution from Dark-Matter-induced photons. We estimate the intensity and features of the angular power spectrum (APS) of this potential Dark Matter (DM) signal, for both decaying and annihilating DM candidates, by constructing template all-sky gamma-ray maps for the emission produced in the galactic halo and its substructures, as well as in extragalactic (sub)halos. The DM distribution is given by state-of-the-art N-body simulations of cosmic structure formation, namely Millennium-II for extragalactic (sub)halos, and Aquarius for the galactic halo and its subhalos. We use a hybrid method of extrapolation to account for (sub)structures that are below the resolution limit of the simulations, allowing us to estimate the total emission all the way down to the minimal self-bound halo mass. We describe in detail the features appearing in the APS of our template maps and we estimate the effect of various uncertainties such as the value of the minimal halo mass, the fraction of substructures hosted in a halo and the shape of the DM density profile. Our results indicate that the fluctuation APS of the DM-induced emission is of the same order as the Fermi-LAT APS, suggesting that one can constrain this hypothetical emission from the comparison with the measured anisotropy. We also quantify the uncertainties affecting our results, finding "theoretical error bands" spanning more than two orders of magnitude and dominated (for a given particle physics model) by our lack of knowledge of the abundance of low-mass (sub)halos.

We present N-body simulations of a new class of self-interacting dark matter models, which do not violate any astrophysical constraints due to a non-power-law velocity dependence of the transfer cross section which is motivated by a Yukawa-like new gauge boson interaction. Specifically, we focus on the formation of a Milky Way-like dark matter halo taken from the Aquarius project and re-simulate it for a couple of representative cases in the allowed parameter space of this new model. We find that for these cases, the main halo only develops a small core (~1 kpc) followed by a density profile identical to that of the standard cold dark matter scenario outside of that radius. Neither the subhalo mass function nor the radial number density of subhaloes are altered in these models but there is a significant change in the inner density structure of subhaloes resulting in the formation of a large density core. As a consequence, the inner circular velocity profiles of the most massive subhaloes differ significantly from the cold dark matter predictions and we demonstrate that they are compatible with the observational data of the brightest Milky Way dSphs in such a velocity-dependent self-interacting dark matter scenario. Specifically, and contrary to the cold dark matter case, there are no subhaloes that are more concentrated than what is inferred from the kinematics of the Milky Way dSphs. We conclude that these models offer an interesting alternative to the cold dark matter model that can reduce the recently reported tension between the brightest Milky Way satellites and the dense subhaloes found in cold dark matter simulations.

We present the first hydrodynamical simulations of structure formation using the new moving mesh code AREPO and compare the results with GADGET simulations based on a traditional smoothed particle hydrodynamics (SPH) technique. The two codes share the same Tree-PM gravity solver and include identical sub-resolution physics, but employ different methods to solve the equations of hydrodynamics. This allows us to assess the impact of hydro-solver uncertainties on the results of cosmological studies of galaxy formation. We focus on predictions for global baryon statistics, such as the cosmic star formation rate density, after we introduce our simulation suite and numerical methods. Properties of individual galaxies and haloes are examined by Keres et al. (2011), while a third paper by Sijacki et al. (2011) uses idealised simulations to analyse the differences between the hydrodynamical schemes. We find that the global baryon statistics differ significantly between the two simulation approaches. AREPO shows higher star formation rates at late times, lower mean temperatures, and different gas mass fractions in characteristic phases of the intergalactic medium, in particular a reduced amount of hot gas. Although both codes use the same implementation of cooling, more gas cools out of haloes in AREPO compared with GADGET towards low redshifts. We show that this is caused by a higher heating rate with SPH in the outer parts of haloes, owing to viscous dissipation of SPH's inherent sonic velocity noise and SPH's efficient damping of subsonic turbulence injected in the halo infall region, and because of a higher efficiency of gas stripping in AREPO. As a result of such differences, AREPO leads also to more disk-like morphologies compared to GADGET. Our results indicate that inaccuracies in hydrodynamic solvers can lead to comparatively large systematic differences.

We present a detailed comparison between the well-known SPH code GADGET and the new moving-mesh code AREPO on a number of hydrodynamical test problems. Through a variety of numerical experiments we establish a clear link between test problems and systematic numerical effects seen in cosmological simulations of galaxy formation. Our tests demonstrate deficiencies of the SPH method in several sectors. These accuracy problems not only manifest themselves in idealized hydrodynamical tests, but also propagate to more realistic simulation setups of galaxy formation, ultimately affecting gas properties in the full cosmological framework, as highlighted in papers by Vogelsberger et al. (2011) and Keres et al. (2011). We find that an inadequate treatment of fluid instabilities in GADGET suppresses entropy generation by mixing, underestimates vorticity generation in curved shocks and prevents efficient gas stripping from infalling substructures. In idealized tests of inside-out disk formation, the convergence rate of gas disk sizes is much slower in GADGET due to spurious angular momentum transport. In simulations where we follow the interaction between a forming central disk and orbiting substructures in a halo, the final disk morphology is strikingly different. In AREPO, gas from infalling substructures is readily depleted and incorporated into the host halo atmosphere, facilitating the formation of an extended central disk. Conversely, gaseous sub-clumps are more coherent in GADGET simulations, morphologically transforming the disk as they impact it. The numerical artefacts of the SPH solver are particularly severe for poorly resolved flows, and thus inevitably affect cosmological simulations due to their hierarchical nature. Our numerical experiments clearly demonstrate that AREPO delivers a physically more reliable solution.

Numerous formulations of finite volume schemes for the Euler and Navier-Stokes equations exist, but in the majority of cases they have been developed for structured and stationary meshes. In many applications, more flexible mesh geometries that can dynamically adjust to the problem at hand and move with the flow in a (quasi) Lagrangian fashion would, however, be highly desirable, as this can allow a significant reduction of advection errors and an accurate realization of curved and moving boundary conditions. Here we describe a novel formulation of viscous continuum hydrodynamics that solves the equations of motion on a Voronoi mesh created by a set of mesh-generating points. The points can move in an arbitrary manner, but the most natural motion is that given by the fluid velocity itself, such that the mesh dynamically adjusts to the flow. Owing to the mathematical properties of the Voronoi tessellation, pathological mesh-twisting effects are avoided. Our implementation considers the full Navier-Stokes equations and has been realized in the AREPO code both in 2D and 3D. We propose a new approach to compute accurate viscous fluxes for a dynamic Voronoi mesh, and use this to formulate a finite volume solver of the Navier-Stokes equations. Through a number of test problems, including circular Couette flow and flow past a cylindrical obstacle, we show that our new scheme combines good accuracy with geometric flexibility, and hence promises to be competitive with other highly refined Eulerian methods. This will in particular allow astrophysical applications of the AREPO code where physical viscosity is important, such as in the hot plasma in galaxy clusters, or for viscous accretion disk models.

We combine the high-resolution Aquarius simulations with three-dimensional models of reionization based on the initial density field of the Aquarius parent simulation, Millennium-II, to study the impact of patchy reionization on the faint satellite population of Milky Way halos. Because the Aquarius suite consists of zoom-in simulations of halos in the Millennium-II volume, we follow the formation of substructure and the growth of reionization bubbles due to the larger environment simultaneously, and thereby determine the reionization redshifts of satellite candidates. We do this for four different reionization models, and also compare results to instantaneous reionization. Using a simple procedure for selecting satellites and assigning luminosities in the simulations, we compare the resulting satellite populations. We find that the overall number of satellites depends sensitively on the reionization model, with a factor of 3-4 variation between the four models for a given host halo, although the difference is entirely in the population of faint satellites (M_V > -10). In addition, we find that for a given reionization model the total number of satellites differs by 10%-20% between the patchy and homogeneous scenarios, provided that the redshift is chosen appropriately for the instantaneous case. However, the halo-halo scatter from the six Aquarius halos is large, up to a factor of 2-3, and so is comparable to the difference between reionization scenarios. In order to use the population of faint dwarf galaxies around the Milky Way as a probe of the local reionization history, then, it is necessary to first better understand the general distribution of substructure around Milky Way-mass halos.

For the first time, the Fermi-LAT measured the angular power spectrum (APS) of anisotropies in the diffuse gamma-ray background. The data is found to be broadly compatible with a model with contributions from the point sources in the 1-year catalog, the galactic diffuse background, and the extragalactic isotropic emission; however deviations are present at both large and small angular scales. In this study, we complement the model with a contribution from Dark Matter (DM) whose distribution is modeled exploiting the results of the most recent N-body simulations, considering the contribution of extragalactic halos and subhalos (from Millennium-II) and of galactic substructures (from Aquarius). With the use of the Fermi Science Tools, these simulations serve as templates to produce mock gamma-ray count maps for DM gamma-ray emission, both in the case of an annihilating and a decaying DM candidate. The APS will then be computed and compared with the Fermi-LAT results to derive constraints on the DM particle physics properties. The possible systematic due to an imperfect model of the galactic foreground is also studied and taken into account properly. The present paper reports on the status of the project.

We use the high resolution cosmological N-body simulations from the Aquarius project to investigate in detail the mechanisms that determine the shape of Milky Way-type dark matter haloes. We find that, when measured at the instantaneous virial radius, the shape of individual haloes changes with time, evolving from a typically prolate configuration at early stages to a more triaxial/oblate geometry at the present day. This evolution in halo shape correlates well with the distribution of the infalling material: prolate configurations arise when haloes are fed through narrow filaments, which characterizes the early epochs of halo assembly, whereas triaxial/oblate configurations result as the accretion turns more isotropic at later times. Interestingly, at redshift z=0, clear imprints of the past history of each halo are recorded in their shapes at different radii, which also exhibit a variation from prolate in the inner regions to triaxial/oblate in the outskirts. Provided that the Aquarius haloes are fair representatives of Milky Way-like 10^12 Msun objects, we conclude that the shape of such dark matter haloes is a complex, time-dependent property, with each radial shell retaining memory of the conditions at the time of collapse.

(Abridged) The extragalactic background light (EBL) observed at multiple wavelengths is a promising tool to probe the nature of dark matter since it might contain a significant contribution from gamma-rays produced promptly by dark matter annihilation. Additionally, the electrons and positrons produced in the annihilation give energy to the CMB photons to populate the EBL with X-rays and gamma-rays. We here create full-sky maps of the radiation from both of these contributions using the high-resolution Millennium-II simulation. We use upper limits on the contributions of unknown sources to the EBL to constrain the intrinsic properties of dark matter using a model-independent approach that can be employed as a template to test different particle physics models (including those with a Sommerfeld enhancement). These upper limits are based on observations spanning eight orders of magnitude in energy (from soft X-rays measured by CHANDRA to gamma-rays measured by Fermi), and on expectations for the contributions from blazars and star forming galaxies. To exemplify this approach, we analyze a set of benchmark Sommerfeld-enhanced models that give the correct dark matter abundance, satisfy CMB constraints, and fit the cosmic ray spectra measured by PAMELA and Fermi without any contribution from local subhalos. We find that these models are in conflict with the EBL constraints unless the contribution of unresolved subhalos is small and the annihilation signal dominates the EBL. We conclude that provided the collisionless cold dark matter paradigm is accurate, even for conservative estimates of the contribution from unresolved subhalos and astrophysical backgrounds, the EBL is at least as sensitive a probe of these types of scenarios as the CMB. Our results disfavor an explanation of the positron excess measured by PAMELA based only on dark matter annihilation in the smooth Galactic halo.

Anomalies in direct and indirect detection have motivated models of dark matter consisting of a multiplet of nearly-degenerate states, coupled by a new GeV-scale interaction. We perform a careful analysis of the thermal freezeout of dark matter annihilation in such a scenario. We compute the range of "boost factors" arising from Sommerfeld enhancement in the local halo for models which produce the correct relic density, and show the effect of including constraints on the saturated enhancement from the cosmic microwave background (CMB). We find that boost factors from Sommerfeld enhancement of up to ~800 are possible in the local halo. When the CMB bounds on the saturated enhancement are applied, the maximal boost factor is reduced to ~400 for 1-2 TeV dark matter and sub-GeV force carriers, but remains large enough to explain the observed Fermi and PAMELA electronic signals. We describe regions in the DM mass-boost factor plane where the cosmic ray data is well fit for a range of final states, and show that Sommerfeld enhancement alone is enough to provide the large annihilation cross sections required to fit the data, although for light mediator masses (less than ~200 MeV) there is tension with the CMB constraints in the absence of astrophysical boost factors from substructure. Additionally, we consider the circumstances under which WIMPonium formation is relevant and find for heavy WIMPs (greater than ~2 TeV) and soft-spectrum annihilation channels it can be an important consideration; we find regions with dark matter mass greater than 2.8 TeV that are consistent with the CMB bounds and have ~600-700 present-day boost factors.

We present the first and so far the only simulations to follow the fine-grained phase-space structure of galaxy haloes formed from generic LCDM initial conditions. We integrate the geodesic deviation equation in tandem with the N-body equations of motion, demonstrating that this can produce numerically converged results for the properties of fine-grained phase-space streams and their associated caustics, even in the inner regions of haloes. Our effective resolution for such structures is many orders of magnitude better than achieved by conventional techniques on even the largest simulations. We apply these methods to the six Milky Way-mass haloes of the Aquarius Project. At 8 kpc from halo centre a typical point intersects about 10^14 streams with a very broad range of individual densities; the ~10^6 most massive streams contribute about half of the local dark matter density. As a result, the velocity distribution of dark matter particles should be very smooth with the most massive fine-grained stream contributing about 0.1% of the total signal. Dark matter particles at this radius have typically passed 200 caustics since the Big Bang. The peak densities on present-day caustics in the inner halo almost all lie well below the mean local dark matter density. As a result caustics provide a negligible boost (<0.1%) to the predicted local dark matter annihilation rate. The effective boost is larger in the outer halo but never exceeds about 10%. Thus fine-grained streams and their associated caustics have no effect on the detectability of dark matter, either directly in Earth-bound laboratories, or indirectly through annihilation radiation, with the exception that resonant cavity experiments searching for axions may see the most massive local fine-grained streams because of their extreme localisation in energy/momentum space. (abridged)

We analyse the coarse-grained phase-space structure of the six Galaxy-scale dark matter haloes of the Aquarius Project using a state-of-the-art 6D substructure finder. Within r_50, we find that about 35% of the mass is in identifiable substructures, predominantly tidal streams, but including about 14% in self-bound subhaloes. The slope of the differential substructure mass function is close to -2, which should be compared to around -1.9 for the population of self-bound subhaloes. Near r_50 about 60% of the mass is in substructures, with about 30% in self-bound subhaloes. The inner 35 kpc of the highest resolution simulation has only 0.5% of its mass in self-bound subhaloes, but 3.3% in detected substructure, again primarily tidal streams. The densest tidal streams near the solar position have a 3-D mass density about 1% of the local mean, and populate the high velocity tail of the velocity distribution.

When a spinning system experiences a transient gravitational encounter with an external perturber, a quasi-resonance occurs if the spin frequency of the victim matches the peak orbital frequency of the perturber. Such encounters are responsible for the formation of long tails and bridges of stars during galaxy collisions. For high-speed encounters, the resulting velocity perturbations can be described within the impulse approximation. The traditional impulse approximation, however, does not distinguish between prograde and retrograde encounters, and therefore completely misses the resonant response. Here, using perturbation theory, we compute the effects of quasi-resonant phenomena on stars orbiting within a disk. Explicit expressions are derived for the velocity and energy change to the stars induced by tidal forces from an external gravitational perturber passing either on a straight line or parabolic orbit. Comparisons with numerical restricted three-body calculations illustrate the applicability of our analysis.

We use the Aquarius simulation series to study the imprint of assembly history on the structure of Galaxy-mass cold dark matter halos. Our results confirm earlier work regarding the influence of mergers on the mass density profile and the inside-out growth of halos. The inner regions that contain the visible galaxies are stable since early times and are significantly affected only by major mergers. Particles accreted diffusely or in minor mergers are found predominantly in the outskirts of halos. Our analysis reveals trends that run counter to current perceptions of hierarchical halo assembly. For example, major mergers (i.e. those with progenitor mass ratios greater than 1:10) contribute little to the total mass growth of a halo, on average less than 20 per cent for our six Aquarius halos. The bulk is contributed roughly equally by minor mergers and by "diffuse" material which is not resolved into individual objects. This is consistent with modeling based on excursion-set theory which suggests that about half of this diffuse material should not be part of a halo of any scale. Interestingly, the simulations themselves suggest that a significantly fraction is not truly diffuse, since it was ejected from earlier halos by mergers prior to their joining the main system. The Aquarius simulations resolve halos to much lower mass scales than are expected to retain gas or form stars. These results thus confirm that most of the baryons from which visible galaxies form are accreted diffusely, rather than through mergers, and they suggest that only relatively rare major mergers will affect galaxy structure at later times.

We carry out fully 3-dimensional simulations of evolution from self-similar, spherically symmetric linear perturbations of a Cold Dark Matter dominated Einstein-de Sitter universe. As a result of the radial orbit instability, the haloes which grow from such initial conditions are triaxial with major-to-minor axis ratios of order 3:1. They nevertheless grow approximately self-similarly in time. In all cases they have power-law density profiles and near-constant velocity anisotropy in their inner regions. Both the power-law index and the value of the velocity anisotropy depend on the similarity index of the initial conditions, the former as expected from simple scaling arguments. Halo structure is thus not "universal" but remembers the initial conditions. On larger scales the density and anisotropy profiles show two characteristic scales, corresponding to particles at first pericentre and at first apocentre after infall. They are well approximated by the NFW model only for one value of the similarity index. In contrast, at all radii within the outer caustic the pseudo phase-space density can be fit by a single power law with an index which depends only very weakly on the similarity index of the initial conditions. This behaviour is very similar to that found for haloes formed from LCDM initial conditions and so can be considered approximately universal.

We calculate how the relic density of dark matter particles is altered when their annihilation is enhanced by the Sommerfeld mechanism due to a Yukawa interaction between the annihilating particles. Maintaining a dark matter abundance consistent with current observational bounds requires the normalization of the s-wave annihilation cross section to be decreased compared to a model without enhancement. The level of suppression depends on the specific parameters of the particle model, with the kinetic decoupling temperature having the most effect. We find that the cross section can be reduced by as much as an order of magnitude for extreme cases. We also compute the mu-type distortion of the CMB energy spectrum caused by energy injection from such Sommerfeld-enhanced annihilation. Our results indicate that in the vicinity of resonances, associated with bound states, distortions can be large enough to be excluded by the upper limit |mu|<9.0x10^(-5) found by the COBE/FIRAS experiment.

We use N-body simulations to investigate the radial dependence of the density and velocity dispersion in cold dark matter (CDM) halos. In particular, we explore how closely Q rho/sigma^3, a surrogate measure of the phase-space density, follows a power-law in radius. Our study extends earlier work by considering, in addition to spherically-averaged profiles, local Q-estimates for individual particles, Q_i; profiles based on the ellipsoidal radius dictated by the triaxial structure of the halo, Q_i(r'); and by carefully removing substructures in order to focus on the profile of the smooth halo, Q^s. The resulting Q_i^s(r') profiles follow closely a power law near the center, but show a clear upturn from this trend near the virial radius, r_{200}. The location and magnitude of the deviations are in excellent agreement with the predictions from Bertschinger's spherical secondary-infall similarity solution. In this model, Q \propto r^{-1.875} in the inner, virialized regions, but departures from a power-law occur near r_{200} because of the proximity of this radius to the location of the first shell crossing - the shock radius in the case of a collisional fluid. Particles there have not yet fully virialized, and so Q departs from the inner power-law profile. Our results imply that the power-law nature of $Q$ profiles only applies to the inner regions and cannot be used to predict accurately the structure of CDM halos beyond their characteristic scale radius.

We simulate the growth of isolated dark matter haloes from self-similar and spherically symmetric initial conditions. Our N-body code integrates the geodesic deviation equation in order to track the streams and caustics associated with individual simulation particles. The radial orbit instability causes our haloes to develop major-to-minor axis ratios approaching 10 to 1 in their inner regions. They grow similarly in time and have similar density profiles to the spherical similarity solution, but their detailed structure is very different. The higher dimensionality of the orbits causes their stream and caustic densities to drop much more rapidly than in the similarity solution. This results in a corresponding increase in the number of streams at each point. At 1% of the turnaround radius (corresponding roughly to the Sun's position in the Milky Way) we find of order 10^6 streams in our simulations, as compared to 10^2 in the similarity solution. The number of caustics in the inner halo increases by a factor of several, because a typical orbit has six turning points rather than one, but caustic densities drop by a much larger factor. This reduces the caustic contribution to the annihilation radiation. For the region between 1% and 50% of the turnaround radius, this is 4% of the total in our simulated haloes, as compared to 6.5% in the similarity solution. Caustics contribute much less at smaller radii. These numbers assume a 100 GeV c^-2 neutralino with present-day velocity dispersion 0.03 cm s^-1, but reducing the dispersion by ten orders of magnitude only doubles the caustic luminosity. We conclude that caustics will be unobservable in the inner parts of haloes. Only the outermost caustic might potentially be detectable.

We study the mass, velocity dispersion, and anisotropy profiles of $\Lambda$CDM halos using a suite of N-body simulations of unprecedented numerical resolution (the {\it Aquarius Project}). Our analysis confirms a number of results claimed by earlier work, and clarifies a few issues where conflicting claims may be found in the recent literature. The spherically-averaged density profile becomes progressively shallower inwards and, at the innermost resolved radius, the logarithmic slope is $\gamma \equiv -$d$\ln\rho/$d$\ln r \simlt 1$. Asymptotic inner slopes as steep as the recently claimed $\rho \propto r^{-1.2}$ are clearly ruled out. The radial dependence of $\gamma$ is well approximated by a power-law, $\gamma \propto r^{\alpha}$ (the Einasto profile). The shape parameter, $\alpha$, varies slightly but significantly from halo to halo, implying that the mass profiles of $\Lambda$CDM halos are not strictly universal: different halos cannot, in general, be rescaled to look identical. Departures from similarity are also seen in velocity dispersion profiles and correlate with those in density profiles so as to preserve a power-law form for the spherically averaged pseudo-phase-space density, $\rho/\sigma^3\propto r^{-1.875}$. Our conclusions are reliable down to radii below 0.4% of the virial radius, providing well-defined predictions for halo structure when baryonic effects are neglected, and thus an instructive theoretical template against which the modifications induced by the baryonic components of real galaxies can be judged.

We study predictions for dark matter phase-space structure near the Sun based on high-resolution simulations of six galaxy halos taken from the Aquarius Project. The local DM density distribution is predicted to be remarkably smooth; the density at the Sun differs from the mean over a best-fit ellipsoidal equidensity contour by less than 15% at the 99.9% confidence level. The local velocity distribution is also very smooth, but it differs systematically from a (multivariate) Gaussian distribution. This is not due to the presence of individual clumps or streams, but to broad features in the velocity modulus and energy distributions that are stable both in space and time and reflect the detailed assembly history of each halo. These features have a significant impact on the signals predicted for WIMP and axion searches. For example, WIMP recoil rates can deviate by ~10% from those expected from the best-fit multivariate Gaussian models. The axion spectra in our simulations typically peak at lower frequencies than in the case of multivariate Gaussian velocity distributions. Also in this case, the spectra show significant imprints of the formation of the halo. This implies that once direct DM detection has become routine, features in the detector signal will allow us to study the dark matter assembly history of the Milky Way. A new field, "dark matter astronomy", will then emerge.

Caustics are a generic feature of the nonlinear growth of structure in the dark matter distribution. If the dark matter were absolutely cold, its mass density would diverge at caustics, and the integrated annihilation probability would also diverge for individual particles participating in them. For realistic dark matter candidates, this behaviour is regularised by small but non-zero initial thermal velocities. We present a mathematical treatment of evolution from Hot, Warm or Cold Dark Matter initial conditions which can be directly implemented in cosmological N-body codes. It allows the identification of caustics and the estimation of their annihilation radiation in fully general simulations of structure formation.

We have performed the largest ever particle simulation of a Milky Way-sized dark matter halo, and present the most comprehensive convergence study for an individual dark matter halo carried out thus far. We have also simulated a sample of 6 ultra-highly resolved Milky-way sized halos, allowing us to estimate the halo-to-halo scatter in substructure statistics. In our largest simulation, we resolve nearly 300,000 gravitationally bound subhalos within the virialized region of the halo. Simulations of the same object differing in mass resolution by factors up to 1800 accurately reproduce the largest subhalos with the same mass, maximum circular velocity and position, and yield good convergence for the abundance and internal properties of dark matter substructures. We detect up to four generations of subhalos within subhalos, but contrary to recent claims, we find less substructure in subhalos than in the main halo when regions of equal mean overdensity are compared. The overall substructure mass fraction is much lower in subhalos than in the main halo. Extrapolating the main halo's subhalo mass spectrum down to an Earth mass, we predict the mass fraction in substructure to be well below 3% within 100 kpc, and to be below 0.1% within the Solar Circle. The inner density profiles of subhalos show no sign of converging to a fixed asymptotic slope and are well fit by gently curving profiles of Einasto form. The mean concentrations of isolated halos are accurately described by the fitting formula of Neto et al. down to maximum circular velocities of 1.5 km/s, an extrapolation over some 5 orders of magnitude in mass. However, at equal maximum circular velocity, subhalos are more concentrated than field halos, with a characteristic density that is typically ~2.6 times larger and increases towards the halo centre.

Dark matter is the dominant form of matter in the universe, but its nature is unknown. It is plausibly an elementary particle, perhaps the lightest supersymmetric partner of known particle species. In this case, annihilation of dark matter in the halo of the Milky Way should produce gamma-rays at a level which may soon be observable. Previous work has argued that the annihilation signal will be dominated by emission from very small clumps (perhaps smaller even than the Earth) which would be most easily detected where they cluster together in the dark matter halos of dwarf satellite galaxies. Here we show, using the largest ever simulation of the formation of a galactic halo, that such small-scale structure will, in fact, have a negligible impact on dark matter detectability. Rather, the dominant and likely most easily detectable signal will be produced by diffuse dark matter in the main halo of the Milky Way. If the main halo is strongly detected, then small dark matter clumps should also be visible, but may well contain no stars, thereby confirming a key prediction of the Cold Dark Matter (CDM) model.

We present a new and completely general technique for calculating the fine-grained phase-space structure of dark matter throughout the Galactic halo. Our goal is to understand this structure on the scales relevant for direct and indirect detection experiments. Our method is based on evaluating the geodesic deviation equation along the trajectories of individual DM particles. It requires no assumptions about the symmetry or stationarity of the halo formation process. In this paper we study general static potentials which exhibit more complex behaviour than the separable potentials studied previously. For ellipsoidal logarithmic potentials with a core, phase mixing is sensitive to the resonance structure, as indicated by the number of independent orbital frequencies. Regions of chaotic mixing can be identified by the very rapid decrease in the real space density of the associated dark matter streams. We also study the evolution of stream density in ellipsoidal NFW halos with radially varying isopotential shape, showing that if such a model is applied to the Galactic halo, at least $10^5$ streams are expected near the Sun. The most novel aspect of our approach is that general non-static systems can be studied through implementation in a cosmological N-body code. Such an implementation allows a robust and accurate evaluation of the enhancements in annihilation radiation due to fine-scale structure such as caustics. We embed the scheme in the current state-of-the-art code GADGET-3 and present tests which demonstrate that N-body discreteness effects can be kept under control in realistic configurations.