We analyse a set of cosmological magneto-hydrodynamic simulations of the formation of Milky Way-mass galaxies identified to have a prominent radially anisotropic stellar halo component similar to the so-called "Gaia Sausage" found in the Gaia data. We examine the effects of the progenitor of the Sausage (the Gaia-Enceladus-Sausage, GES) on the formation of major galactic components analogous to the Galactic thick disc and inner stellar halo. We find that the GES merger is likely to have been gas-rich and contribute 10-50$\%$ of gas to a merger-induced centrally concentrated starburst that results in the rapid formation of a compact, rotationally supported thick disc that occupies the typical chemical thick disc region of chemical abundance space. We find evidence that gas-rich mergers heated the proto-disc of the Galaxy, scattering stars onto less-circular orbits such that their rotation velocity and metallicity positively correlate, thus contributing an additional component that connects the Galactic thick disc to the inner stellar halo. We demonstrate that the level of kinematic heating of the proto-galaxy correlates with the kinematic state of the population before the merger, the progenitor mass and orbital eccentricity of the merger. Furthermore, we show that the mass and time of the merger can be accurately inferred from local stars on counter-rotating orbits.

We study the properties of the dark matter component of the radially anisotropic stellar population recently identified in the Gaia data, using magneto-hydrodynamical simulations of Milky Way-like halos from the Auriga project. We identify 10 simulated galaxies that approximately match the rotation curve and stellar mass of the Milky Way. Four of these have an anisotropic stellar population reminiscent of the Gaia structure. We find an anti-correlation between the dark matter mass fraction of this population in the Solar neighbourhood and its orbital anisotropy. We estimate the local dark matter density and velocity distribution for halos with and without the anisotropic stellar population, and use them to simulate the signals expected in future xenon and germanium direct detection experiments. We find that a generalized Maxwellian distribution fits the dark matter halo integrals of the Milky Way-like halos containing the radially anisotropic stellar population. For dark matter particle masses below approximately 10 GeV, direct detection exclusion limits for the simulated halos with the anisotropic stellar population show a mild shift towards smaller masses compared to the commonly adopted Standard Halo Model.

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.

We investigate whether the galaxy and star formation model used for the Auriga simulations can produce a realistic globular cluster (GC) population. We compare statistics of GC candidate star particles in the Auriga haloes with catalogues of the Milky Way (MW) and Andromeda (M31) GC populations. We find that the Auriga simulations do produce sufficient stellar mass for GC candidates at radii and metallicities that are typical for the MW GC system (GCS). We also find varying mass-ratios of the simulated GC candidates relative to the observed mass in the MW and M31 GC systems for different bins of galactocentric radius-metallicity (r$_{\text{gal}}$ -[Fe/H]). Overall, the Auriga simulations produce GC candidates with higher metallicities than the MW and M31 GCS and they are found at larger radii than observed. The Auriga simulations would require bound cluster formation efficiencies higher than ten percent for the metal-poor GC candidates, and those within the Solar radius should experience negligible destruction rates to be consistent with observations. GC candidates in the outer halo, on the other hand, should either have low formation efficiencies, or experience high mass loss for the Auriga simulations to produce a GCS that is consistent with that of the MW or M31. Finally, the scatter in the metallicity as well as in the radial distribution between different Auriga runs is considerably smaller than the differences between that of the MW and M31 GCSs. The Auriga model is unlikely to give rise to a GCS that can be consistent with both 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.

We explore the chemodynamical properties of a sample of barred galaxies in the Auriga magneto-hydrodynamical cosmological zoom-in simulations, which form boxy/peanut (b/p) bulges, and compare these to the Milky Way (MW). We show that the Auriga galaxies which best reproduce the chemodynamical properties of stellar populations in the MW bulge have quiescent merger histories since redshift $z\sim3.5$: their last major merger occurs at $t_{\rm lookback}>12\,\rm Gyrs$, while subsequent mergers have a stellar mass ratio of $\leq$1:20, suggesting an upper limit of a few percent for the mass ratio of the recently proposed Gaia Sausage/Enceladus merger. These Auriga MW-analogues have a negligible fraction of ex-situ stars in the b/p region ($<1\%$), with flattened, thick disc-like metal-poor stellar populations. The average fraction of ex-situ stars in the central regions of all Auriga galaxies with b/p's is 3% -- significantly lower than in those which do not host a b/p or a bar. While the central regions of these barred galaxies contain the oldest populations, they also have stars younger than 5Gyrs (>30%) and exhibit X-shaped age and abundance distributions. Examining the discs in our sample, we find that in some cases a star-forming ring forms around the bar, which alters the metallicity of the inner regions of the galaxy. Further out in the disc, bar-induced resonances lead to metal-rich ridges in the $V_{\phi}-r$ plane -- the longest of which is due to the Outer Lindblad Resonance. Our results suggest the Milky Way has an uncommonly quiet merger history, which leads to an essentially in-situ bulge, and highlight the significant effects the bar can have on the surrounding disc.

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 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 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 use cosmological, magnetohydrodynamical simulations of Milky Way-mass galaxies from the Auriga project to study their enrichment with rapid neutron capture (r-process) elements. We implement a variety of enrichment models from both binary neutron star mergers and rare core-collapse supernovae. We focus on the abundances of (extremely) metal-poor stars, most of which were formed during the first ~Gyr of the Universe in external galaxies and later accreted onto the main galaxy. We find that the majority of metal-poor stars are r-process enriched in all our enrichment models. Neutron star merger models result in a median r-process abundance ratio which increases with metallicity, whereas the median trend in rare core-collapse supernova models is approximately flat. The scatter in r-process abundance increases for models with longer delay times or lower rates of r-process producing events. Our results are nearly perfectly converged, in part due to the mixing of gas between mesh cells in the simulations. Additionally, different Milky Way-mass galaxies show only small variation in their respective r-process abundance ratios. Current (sparse and potentially biased) observations of metal-poor stars in the Milky Way seem to prefer rare core-collapse supernovae over neutron star mergers as the dominant source of r-process elements at low metallicity, but we discuss possible caveats to our models. Dwarf galaxies which experience a single r-process event early in their history show highly enhanced r-process abundances at low metallicity, which is seen both in observations and in our simulations. We also find that the elements produced in a single event are mixed with ~10^8 Msun of gas relatively quickly, distributing the r-process elements over a large region.

We use magneto-hydrodynamical simulations of Milky Way-mass haloes from the Auriga project to examine the properties of surviving and destroyed dwarf galaxies that are accreted by these haloes over cosmic time. We show that the combined luminosity function of surviving and destroyed dwarfs at infall is similar in the various Auriga haloes, and is dominated by the destroyed dwarfs. There is, however, a strong dependence on infall time: destroyed dwarfs have typically early infall times, $t_{infall}<6$ Gyr, whereas the majority of dwarfs accreted at $t_{infall}>10$ Gyr have survived to the present day. Because of their late infall the surviving satellites today had higher metallicites at infall than their destroyed counterparts of similar infall mass; the difference is even more pronounced for the present-day metallicites of satellites, many of which continue to form stars after infall. In agreement with previous work, we find that a small number of relatively massive destroyed dwarf galaxies dominate the mass of the stellar haloes. However, there is a significant radial dependence: while 90 per cent of the mass in the inner regions ($<\,20\,$kpc) is contributed, on average, by only 3 massive progenitors, the outer regions ($>\,100\,$kpc) typically have $\sim8$ main progenitors of relatively lower mass. Finally, we show that a few massive progenitors dominate the metallicity distribution of accreted stars, even at the metal poor end. Contrary to common assumptions in the literature, dwarf galaxies of mass $M_{*}<10^7 \, M_{\odot}$ make up less than 10 per cent of the accreted, metal poor stars ([Fe/H] $<\,-3$) in the inner $50\,$kpc.

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 study the orbital phase-space of dark matter (DM) halos in the AURIGA suite of cosmological hydrodynamics simulations of Milky Way analogues. We characterise halos by their spherical action distribution, a function of the specific angular momentum, and the radial action, of the DM particles. By comparing DM-only and hydrodynamical simulations of the same halos, we investigate the contraction of DM halos caused by the accumulation of baryons at the centre. We find a small systematic suppression of the radial action in the DM halos of the hydrodynamical simulations, suggesting that the commonly used adiabatic contraction approximation can result in an underestimate of the density by ~ 8%. We apply an iterative algorithm to contract the AURIGA DM halos given a baryon density profile and halo mass, recovering the true contracted DM profiles with an accuracy of ~15%, that reflects halo-to-halo variation. Using this algorithm, we infer the total mass profile of the Milky Way's contracted DM halo. We derive updated values for the key astrophysical inputs to DM direct detection experiments: the DM density and velocity distribution in the Solar neighbourhood.

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.

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.

The circumgalactic medium (CGM) is one of the frontiers of galaxy formation and intimately connected to the galaxy via accretion of gas on to the galaxy and gaseous outflows from the galaxy. Here we analyse the magnetic field in the CGM of the Milky Way-like galaxies simulated as part of the \textsc{Auriga} project that constitutes a set of high resolution cosmological magnetohydrodynamical zoom simulations. We show that at high redshift the CGM becomes magnetised via galactic outflows that transport magnetised gas from the disk into the halo. At this time the magnetisation of the CGM closely follows its metal enrichment. We then show that at low redshift an in-situ turbulent dynamo that operates on a timescale of Gigayears further amplifies the magnetic field in the CGM and sets its strength down $z=0$. The magnetic field strength reaches a typical value of $0.1\,\mu G$ at the virial radius at $z=0$ and becomes mostly uniform within the virial radius. Its Faraday rotation signal is in excellent agreement with recent observations. For most of its evolution the magnetic field in the CGM is an unordered small scale field. Only strong coherent outflows at low redshift are able to order the magnetic field in parts of the CGM that are directly displaced by these outflows.

We use the Auriga cosmological simulations of Milky Way (MW)-mass galaxies and their surroundings to study the satellite populations of dwarf galaxies in $\Lambda$CDM. As expected from prior work, the number of satellites above a fixed stellar mass is a strong function of the mass of the primary dwarf. For galaxies as luminous as the Large Magellanic Cloud (LMC), and for halos as massive as expected for the LMC (determined by its rotation speed), the simulations predict about 3 satellites with stellar masses exceeding $M_*>10^5\, M_\odot$. If the LMC is on its first pericentric passage, then these satellites should be near the LMC and should have orbital angular momenta roughly coincident with that of the LMC. We use 3D positions and velocities from the 2nd data release of the Gaia mission to revisit which of the "classical" MW dwarf spheroidals could plausibly be LMC satellites. The new proper motions of the Fornax and Carina dwarf spheroidals place them on orbits closely aligned with the orbital plane of the Magellanic Clouds, hinting at a potential Magellanic association. Together with the Small Magellanic Cloud (SMC), this result raises to $3$ the number of LMC satellites with $M_*>10^5\, M_\odot$, as expected from simulations. This also fills the 12-mag luminosity gap between the SMC and the ultra-faints Hyi1, Car2, Hor1, and Car3, the few ultra-faint satellites confirmed to have orbits consistent with a Magellanic origin.

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.

In this work we analyse the structural and photometric properties of 21 barred simulated galaxies from the Auriga Project. These consist of Milky Way-mass magneto-hydrodynamical simulations in a $\Lambda$CDM cosmological context. In order to compare with observations, we generate synthetic SDSS-like broad-band images from the numerical data at z = 0 with different inclinations (from face-on to edge-on). Ellipse fits are used to determine the bar lengths, and 2D bulge/disc/bar decompositions with galfit are also performed, modelling the bar component with the modified Ferrer profile. We find a wide range of bar sizes and luminosities in the sample, and their structural parameters are in good agreement with the observations. All bulges present low S\'ersic indexes, and are classified as pseudobulges. In regard to the discs, the same breaks in the surface brightness profiles observed in real galaxies are found, and the radii at which these take place are in agreement with the observations. Also, from edge-on unsharp-masked images at z = 0, boxy or peanut-shaped (B/P) structures are clearly identified in the inner part of 4 bars, and also 2 more bars are found in buckling phase. The sizes of the B/P match fairly well with those obtained from observations. We thus conclude that the observed photometric and structural properties of galaxies with bars, which are the main drivers of secular evolution, can be developed in present state-of-the-art $\Lambda$CDM cosmological simulations.

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.

We investigate the formation of ultra-diffuse galaxies (UDGs) using the Auriga high-resolution cosmological magneto-hydrodynamical simulations of Milky Way-sized galaxies. We identify a sample of $92$ UDGs in the simulations that match a wide range of observables such as sizes, central surface brightness, S\'{e}rsic indices, colors, spatial distribution and abundance. Auriga UDGs have dynamical masses similar to normal dwarfs. In the field, the key to their origin is a strong correlation present in low-mass dark matter haloes between galaxy size and halo spin parameter. Field UDGs form in dark matter haloes with larger spins compared to normal dwarfs in the field, in agreement with previous semi-analytical models. Satellite UDGs, on the other hand, have two different origins: $\sim 55\%$ of them formed as field UDGs before they were accreted; the remaining $\sim 45\%$ were normal field dwarfs that subsequently turned into UDGs as a result of tidal interactions.

Using a set of 15 high-resolution magnetohydrodynamic cosmological simulations of Milky Way formation, we investigate the origin of the baryonic material found in stars at redshift zero. We find that roughly half of this material originates from subhalo/satellite systems and half is smoothly accreted from the Inter-Galactic Medium (IGM). About $90 \%$ of all material has been ejected and re-accreted in galactic winds at least once. The vast majority of smoothly accreted gas enters into a galactic fountain that extends to a median galactocentric distance of $\sim 20$ kpc with a median recycling timescale of $\sim 500$ Myr. We demonstrate that, in most cases, galactic fountains acquire angular momentum via mixing of low-angular momentum, wind-recycled gas with high-angular momentum gas in the Circum-Galactic Medium (CGM). Prograde mergers boost this activity by helping to align the disc and CGM rotation axes, whereas retrograde mergers cause the fountain to lose angular momentum. Fountain flows that promote angular momentum growth are conducive to smooth evolution on tracks quasi-parallel to the disc sequence of the stellar mass-specific angular momentum plane, whereas retrograde minor mergers, major mergers and bar-driven secular evolution move galaxies towards the bulge-sequence. Finally, we demonstrate that fountain flows act to flatten and narrow the radial metallicity gradient and metallicity dispersion of disc stars, respectively. Thus, the evolution of galactic fountains depends strongly on the cosmological merger history and is crucial for the chemo-dynamical evolution of Milky Way-sized disc 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.

Reionization leads to large spatial fluctuations in the intergalactic temperature that can persist well after its completion. We study the imprints of such fluctuations on the $z\sim5$ Ly$\alpha$ forest flux power spectrum using a set of radiation-hydrodynamic simulations that model different reionization scenarios. We find that large-scale coherent temperature fluctuations bring $\sim20-60\%$ extra power at $k\sim0.002$ s/km, with the largest enhancements in the models where reionization is extended or ends the latest. On smaller scales ($k\gtrsim0.1$ s/km), we find that temperature fluctuations suppress power by $\lesssim10\%$. We find that the shape of the power spectrum is mostly sensitive to the reionization midpoint rather than temperature fluctuations from reionization's patchiness. However, for all of our models with reionization midpoints of $z\le 8$ ($z\le 12$) the shape differences are $\lesssim20\%$ ($\lesssim40\%$) because of a surprisingly well-matched cancellation between thermal broadening and pressure smoothing that occurs for realistic thermal histories. We also consider fluctuations in the ultraviolet background, finding their impact on the power spectrum to be much smaller than temperature fluctuations at $k\gtrsim0.01$ s/km. Furthermore, we compare our models to power spectrum measurements, finding that none of our models with reionization midpoints of $z<8$ is strongly preferred over another and that all of our models with midpoints of $z\geq8$ are excluded at $2.5\sigma$. Future measurements may be able to distinguish between viable reionization models if they can be performed at lower $k$ or, alternatively, if the error bars on the high-$k$ power can be reduced by a factor of $1.5$.

We analytically derive the observed size-mass relation of galaxies' atomic hydrogen (HI), including limits on its scatter, based on simple assumptions about the structure of HI discs. We trial three generic profiles for HI surface density as a function of radius. Firstly, we assert that HI surface densities saturate at a variable threshold, and otherwise fall off exponentially with radius or, secondly, radius squared. Our third model assumes the total gas surface density is exponential, with the HI fraction at each radius depending on local pressure. These are tested against a compilation of 110 galaxies from the THINGS, LITTLE THINGS, LVHIS, and Bluedisk surveys, whose HI surface density profiles are well resolved. All models fit the observations well and predict consistent size-mass relations. Using an analytical argument, we explain why processes that cause gas disc truncation - such as ram-pressure stripping - scarcely affect the HI size-mass relation. This is tested with the IllustrisTNG(100) cosmological, hydrodynamic simulation and the Dark Sage semi-analytic model of galaxy formation, both of which capture radially resolved disc structure. For galaxies with m_*>10^9 M_solar and m_HI>10^8 M_solar, both simulations predict HI size-mass relations that align with observations, show no difference between central and satellite galaxies, and show only a minor, second-order dependence on host halo mass for satellites. Ultimately, the universally tight HI size-mass relation is mathematically inevitable and robust. Only by completely disrupting the structure of HI discs, e.g. through overly powerful feedback, could a simulation predict the relation poorly.

We explore the predictive power of cosmological, hydrodynamical simulations for stellar phase space substructure and velocity correlations with the Auriga simulations and Aurigaia mock-Gaia catalogues. We show that at the solar circle the Auriga simulations commonly host phase space structures in the stellar component that have constant orbital energies and arise from accreted subhaloes. These structures can persist for a few Gyrs, even after coherent streams in position space have been erased. We also explore velocity two-point correlation functions and find this diagnostic is not deterministic for particular clustering patterns in phase space. Finally, we explore these structure diagnostics with the Aurigaia catalogues and show that the current catalogues have the ability to recover some structures in phase space but careful consideration is required to separate physical structures from numerical structures arising from catalogue generation methods.

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.

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 analyse the orbital kinematics of the Milky Way (MW) satellite system utilizing the latest systemic proper motions for 38 satellites based on data from Gaia Data Release 2. Combining these data with distance and line-of-sight velocity measurements from the literature, we use a likelihood method to model the velocity anisotropy, $\beta$, as a function of Galactocentric distance and compare the MW satellite system with those of simulated MW-mass haloes from the APOSTLE and Auriga simulation suites. The anisotropy profile for the MW satellite system increases from $\beta\sim -2$ at $r\sim20$ kpc to $\beta\sim 0.5$ at $r\sim200$ kpc, indicating that satellites closer to the Galactic centre have tangentially-biased motions while those farther out have radially-biased motions. The motions of satellites around APOSTLE host galaxies are nearly isotropic at all radii, while the $\beta(r)$ profiles for satellite systems in the Auriga suite, whose host galaxies are substantially more massive in baryons than those in APOSTLE, are more consistent with that of the MW satellite system. This shape of the $\beta(r)$ profile may be attributed to the central stellar disc preferentially destroying satellites on radial orbits, or intrinsic processes from the formation of the Milky Way system.

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.

N-body simulations make unambiguous predictions for the abundance of substructures within dark matter halos. However, the inclusion of baryons in the simulations changes the picture because processes associated with the presence of a large galaxy in the halo can destroy subhalos and substantially alter the mass function and velocity distribution of subhalos. We compare the effect of galaxy formation on subhalo populations in two state-of-the-art sets of hydrodynamical CDM simulations of Milky Way mass halos, APOSTLE and AURIGA. We introduce a new method for tracking the orbits of subhalos between simulation snapshots that gives accurate results down to a few kiloparsecs from the centre of the halo. Relative to a dark matter-only simulation, the abundance of subhalos in APOSTLE is reduced by 50% near the centre and by 10% within r200. In AURIGA the corresponding numbers are 80% and 40%. The velocity distributions of subhalos are also affected by the presence of the galaxy, much more so in AURIGA than in APOSTLE . The differences on subhalo properties in the two simulations can be traced back to the mass of the central galaxies, which in AURIGA are typically twice as massive as those in APOSTLE . We show that some of the results from previous studies are inaccurate due to systematic errors in the modelling of subhalo orbits near the centre of halos.

We study the galactic bulges in the Auriga simulations, a suite of thirty cosmological magneto-hydrodynamical zoom-in simulations of late-type galaxies in Milky Way-sized dark matter haloes performed with the moving-mesh code AREPO. We aim to characterize bulge formation mechanisms in this large suite of galaxies simulated at high resolution in a fully cosmological context. The bulges of the Auriga galaxies show a large variety in their shapes,sizes and formation histories. According to observational classification criteria, such as Sersic index and degree of ordered rotation, the majority of the Auriga bulges can be classified as pseudo-bulges, while some of them can be seen as composite bulges with a classical component; however, none can be classified as a classical bulge. Auriga bulges show mostly an in-situ origin, 21 percent of them with a negligible accreted fraction (facc < 0.01). In general,their in-situ component was centrally formed, with 75 percent of the bulges forming most of their stars inside the bulge region at z=0. Part of their in-situ mass growth is rapid and is associated with the effects of mergers, while another part is more secular in origin. In 90 percent of the Auriga bulges, the accreted bulge component originates from less than four satellites.We investigate the relation between the accreted stellar haloes and the bulges of the Auriga simulations. The total bulge mass shows no correlation with the accreted stellar halo mass, as in observations. However, the accreted mass of bulges tends to correlate with their respective accreted stellar halo mass.

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.

The dark matter velocity distribution in the Solar neighbourhood is an important astrophysical input which enters in the predicted event rate of dark matter direct detection experiments. It has been recently suggested that the local dark matter velocity distribution can be inferred from that of old or metal-poor stars in the Milky Way. We investigate this potential relation using six high resolution magneto-hydrodynamical simulations of Milky Way-like galaxies of the Auriga project. We do not find any correlation between the velocity distributions of dark matter and old stars in the Solar neighbourhood. Likewise, there are no strong correlations between the local velocity distributions of dark matter and metal-poor stars selected by applying reasonable cuts on metallicity. In some simulated galaxies, extremely metal-poor stars have a velocity distribution that is statistically consistent with that of the dark matter, but the sample of such stars is so small that we cannot draw any strong conclusions.

We investigate the impact of dynamical streams and substructure on estimates of the local escape speed and total mass of Milky Way-mass galaxies from modelling the high velocity tail of local halo stars. We use a suite of high-resolution, magneto-hydrodynamical cosmological zoom-in simulations, which resolve phase space substructure in local volumes around solar-like positions. We show that phase space structure varies significantly between positions in individual galaxies and across the suite. Substructure populates the high velocity tail unevenly and leads to discrepancies in the mass estimates. We show that a combination of streams, sample noise and truncation of the high velocity tail below the escape speed leads to a distribution of mass estimates with a median that falls below the true value by $\sim 20 \%$, and a spread of a factor of 2 across the suite. Correcting for these biases, we derive a revised value for the Milky Way mass presented in Deason et al. of $1.29 ^{+0.37}_{-0.47} \times 10^{12}$ $\rm M_{\odot}$.

Measurements of the rotation curves of dwarf galaxies are often interpreted as requiring a constant density core at the centre, at odds with the "cuspy" inner profiles predicted by $N$-body simulations of cold dark matter (CDM) haloes. It has been suggested that this conflict could be resolved by fluctuations in the inner gravitational potential caused by the periodic removal of gas following bursts of star formation. Earlier work has suggested that core formation requires a bursty and extended star formation history (SFH). Here we investigate the structure of CDM haloes of dwarf galaxies ($M_{{\rm DM}} \sim 10^9-5\times10^{10}\,{\rm M}_\odot$) formed in the APOSTLE ('A Project of Simulating the Local Environment') and AURIGA cosmological hydrodynamic simulations. Our simulations have comparable or better resolution than others that make cores ($M_{{\rm gas}} \sim 10^4\,{\rm M}_\odot$, gravitational softening $\sim 150$ pc). Yet, we do not find evidence of core formation at {\it any} mass or any correlation between the inner slope of the DM density profile and temporal variations in the SFH. APOSTLE and AURIGA dwarfs display a similar diversity in their cumulative SFHs to available data for Local Group dwarfs. Dwarfs in both simulations are DM-dominated on all resolved scales at all times, likely limiting the ability of gas outflows to alter significantly the central density profiles of their haloes. We conclude that recurrent bursts of star formation are not sufficient to cause the formation of cores, and that other conditions must also be met for baryons to be able to modify the central DM cusp.

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 have carried out a controlled comparison of the structural and kinematic properties of the circumgalactic medium (CGM) around Milky Way mass galaxies in the Illustris and IllustrisTNG simulations. Very striking differences are found. At z=0, gas column density and temperature profiles at large radii ($\sim 100$ kpc) correlate strongly with disk gas mass fraction in Illustris, but not in TNG. The neutral gas at large radii is preferentially aligned in the plane of the disk in TNG, whereas it is much more isotropic in Illustris. The vertical coherence scale of the rotationally supported gas in the CGM is linked to the gas mass fraction of the galaxy in Illustris, but not in TNG. A tracer particle analysis allows us to show how these differences can be understood as a consequence of the different sub-grid models of feedback in the two simulations. A study of spatially matched galaxies in the two simulations shows that in TNG, feedback by supernovae and AGN helps to create an extended smooth reservoir of hot gas at high redshifts, that then cools to form a thin, rotationally-supported disk at later times. In Illustris, AGN dump heat in the form of hot gas bubbles that push diffuse material at large radii out of the halo. The disk is formed by accretion of colder, recycled material, and this results in more vertically extended gas distributions above and below the Galactic plane. We conclude that variations in the structure of gas around Milky Way mass galaxies are a sensitive probe of feedback physics in simulations and are worthy of more observational consideration in future.

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 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 model the fastest moving (v_tot > 300 km/s) local (D < 3 kpc) halo stars using cosmological simulations and 6-dimensional Gaia data. Our approach is to use our knowledge of the assembly history and phase-space distribution of halo stars to constrain the form of the high velocity tail of the stellar halo. Using simple analytical models and cosmological simulations, we find that the shape of the high velocity tail is strongly dependent on the velocity anisotropy and number density profile of the halo stars --- highly eccentric orbits and/or shallow density profiles have more extended high velocity tails. The halo stars in the solar vicinity are known to have a strongly radial velocity anisotropy, and it has recently been shown the origin of these highly eccentric orbits is the early accretion of a massive (M_star ~ 10^9 M_Sun) dwarf satellite. We use this knowledge to construct a prior on the shape of the high velocity tail. Moreover, we use the simulations to define an appropriate outer boundary of 2r_200, beyond which stars can escape. After applying our methodology to the Gaia data, we find a local (r_0=8.3 kpc) escape speed of v_esc(r_0) = 528(+24,-25) km/s. We use our measurement of the escape velocity to estimate the total Milky Way mass, and dark halo concentration: M_200,tot = 1.00(+0.31,-0.24) x 10^12 M_Sun, c_200 = 10.9(+4.4,-3.3). Our estimated mass agrees with recent results in the literature that seem to be converging on a Milky Way mass of M_200,tot ~ 10^12 M_Sun.

We examine the stellar haloes of the Auriga simulations, a suite of thirty cosmological magneto-hydrodynamical high-resolution simulations of Milky Way-mass galaxies performed with the moving-mesh code AREPO. We study halo global properties and radial profiles out to $\sim 150$ kpc for each individual galaxy. The Auriga haloes are diverse in their masses and density profiles; mean metallicity and metallicity gradients; ages; and shapes, reflecting the stochasticity inherent in their accretion and merger histories. A comparison with observations of nearby late-type galaxies shows very good agreement between most observed and simulated halo properties. However, Auriga haloes are typically too massive. We find a connection between population gradients and mass assembly history: galaxies with few significant progenitors have more massive haloes, possess large negative halo metallicity gradients and steeper density profiles. The number of accreted galaxies, either disrupted or under disruption, that contribute 90% of the accreted halo mass ranges from 1 to 14, with a median of 6.5, and their stellar masses span over three orders of magnitude. The observed halo mass--metallicity relation is well reproduced by Auriga and is set by the stellar mass and metallicity of the dominant satellite contributors. This relationship is found not only for the accreted component but also for the total (accreted + in-situ) stellar halo. Our results highlight the potential of observable halo properties to infer the assembly history of galaxies.

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.

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 present and apply a method to infer the mass of the Milky Way (MW) by comparing the dynamics of MW satellites to those of model satellites in the EAGLE cosmological hydrodynamics simulations. A distribution function (DF) for galactic satellites is constructed from EAGLE using specific angular momentum and specific energy, which are scaled so as to be independent of host halo mass. In this 2-dimensional space, the orbital properties of satellite galaxies vary according to the host halo mass. The halo mass can be inferred by calculating the likelihood that the observed satellite population is drawn from this DF. Our method is robustly calibrated on mock EAGLE systems. We validate it by applying it to the completely independent suite of 30 AURIGA high-resolution simulations of MW-like galaxies: the method accurately recovers their true mass and associated uncertainties. We then apply it to ten classical satellites of the MW with 6D phase-space measurements, including updated proper motions from the GAIA satellite. The mass of the MW is estimated to be $M_{200}^{\textnormal{MW}}=1.17_{-0.15}^{+0.21}\times10^{12}M_{\odot}$ (68\% confidence limits). We combine our total mass estimate with recent mass estimates in the inner regions of the Galaxy to infer an inner dark matter (DM) mass fraction $M^\textnormal{DM}(<20~\rm{kpc})/M^\textnormal{DM}_{200}=0.12$ which is typical of ${\sim}10^{12}M_{\odot}$ $\Lambda$CDM haloes in hydrodynamical galaxy formation simulations. Assuming an NFW profile, this is equivalent to a halo concentration of $c_{200}^{\textnormal{MW}}=10.9^{+2.6}_{-2.0}$.

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 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.

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 present a suite of high resolution radiation hydrodynamic simulations of a small patch ($1 \ {\rm kpc}^2$) of the inter-stellar medium (ISM) performed with Arepo-RT, with the aim to quantify the efficacy of various feedback processes like supernovae explosions (SNe), photoheating and radiation pressure in low gas surface density galaxies ($\Sigma_{\rm gas} \simeq 10 \ {\rm M}_\odot \ {\rm pc}^{-2}$). We show that radiation fields decrease the star formation rate and therefore the total stellar mass formed by a factor of $\sim 2$. This increases the gas depletion timescale and brings the simulated Kennicutt-Schmidt relation closer to the observational estimates. Radiation feedback coupled with SNe is more efficient at driving outflows with the mass and energy loading increasing by a factor of $\sim 10$. This increase is mainly driven by the additional entrainment of medium density ($10^{-2} \leq n< 1 \ {\rm cm}^{-3}$), warm ($300 \ {\rm K}\leq T<8000 \ {\rm K}$) material. Therefore including radiation fields tends to launch colder, denser and higher mass and energy loaded outflows. This is because photoheating of the high density gas around a newly formed star over-pressurises the region, causing it to expand. This reduces the ambient density in which the SNe explode by a factor of $10-100$ which in turn increases their momentum output by a factor of $\sim 1.5-2.5$. Finally, we note that in these low gas surface density environments, radiation fields primarily impact the ISM via photoheating and radiation pressure has only a minimal role in regulating star formation.

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.

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.

Using the astrometry from the ESA's Gaia mission, previous works have shown that the Milky Way stellar halo is dominated by metal-rich stars on highly eccentric orbits. To shed light on the nature of this prominent halo component, we have analysed 28 Galaxy analogues in the Auriga suite of cosmological hydrodynamics zoom-in simulations. Some three quarters of the Auriga galaxies contain significant components with high radial velocity anisotropy, beta > 0.6. However, only in one third of the hosts do the high-beta stars contribute significantly to the accreted stellar halo overall, similar to what is observed in the Milky Way. For this particular subset we reveal the origin of the dominant stellar halo component with high metallicity, [Fe/H]~-1, and high orbital anisotropy, beta>0.8, by tracing their stars back to the epoch of accretion. It appears that, typically, these stars come from a single dwarf galaxy with a stellar mass of order of 10^9-10^10 Msol that merged around 6-10 Gyr ago, causing a sharp increase in the halo mass. Our study therefore establishes a firm link between the excess of radially anisotropic stellar debris in the Milky Way halo and an ancient head-on collision between the young Milky Way and a massive dwarf galaxy

Faraday rotation is one of the most widely used observables to infer the strength and configuration of the magnetic field in the ionised gas of the Milky Way and nearby spiral galaxies. Here we compute synthetic Faraday rotation maps at $z=0$ for a set of disk galaxies from the Auriga high-resolution cosmological simulations, for different observer positions within and outside the galaxy. We find that the strength of the Faraday rotation of our simulated galaxies for a hypothetic observer at the solar circle is broadly consistent with the Faraday rotation seen for the Milky Way. The same holds for an observer outside the galaxy and the observed signal of the nearby spiral galaxy M51. However, we also find that the structure and angular power spectra of the synthetic all-sky Faraday rotation maps vary strongly with azimuthal position along the solar circle. We argue that this variation is a result of the structure of the magnetic field of the galaxy that is dominated by an azimuthal magnetic field ordered scales of several kpc, but has radial and vertical magnetic field components that are only ordered on scales of 1-2 kpc. Because the magnetic field strength decreases exponentially with height above the disk, the Faraday rotation for an observer at the solar circle is dominated by the local environment. This represents a severe obstacle for attempts to reconstruct the global magnetic field of the Milky Way from Faraday rotation maps alone without including additional observables.

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.

We present and analyse mock stellar catalogues that match the selection criteria and observables (including uncertainties) of the Gaia satellite data release 2 (DR2). The source are six cosmological high-resolution magneto-hydrodynamic $\Lambda$CDM zoom simulations of the formation of Milky Way analogues from the AURIGA project. Mock data are provided for stars with $V < 16$ mag, and $V < 20$ mag at $|b|>20$ degrees. The mock catalogues are made using two different methods: the public SNAPDRAGONS code, and a method based on that of Lowing et al. that preserves the phase-space distribution of the model stars. These publicly available catalogues contain 5-parameter astrometry, radial velocities, multi-band photometry, stellar parameters, dust extinction values, and uncertainties in all these quantities. In addition, we provide the gravitational potential and information on the origin of each star. By way of demonstration, we apply the mock catalogues to analyses of the young stellar disc and the stellar halo. We show that: i) the young outer stellar disc exhibits a flared distribution that is detectable in the height and vertical velocity distribution of A- and B-dwarf stars up to radii of ~15 kpc; and ii) the spin of the stellar halo out to 100 kpc can be accurately measured with Gaia DR2 RR Lyrae stars. These catalogues are well suited for comparisons with observations and should help to: i) develop and test analysis methods for the Gaia DR2 data; ii) gauge the limitations and biases of the data and iii) interpret the data in the light of theoretical predictions from realistic $ab$ $initio$ simulations of galaxy formation in the $\Lambda$CDM cosmological model.

We present a suite of six fully cosmological, three-dimensional simulations of the collapse of an atomic cooling halo in the early Universe. We use the moving-mesh code arepo with an improved primordial chemistry network to evolve the hydrodynamical and chemical equations. The addition of a strong Lyman-Werner background suppresses molecular hydrogen cooling and permits the gas to evolve nearly isothermally at a temperature of about 8000 K. Strong gravitational torques effectively remove angular momentum and lead to the central collapse of gas, forming a supermassive protostar at the center of the halo. We model the protostar using two methods: sink particles that grow through mergers with other sink particles, and a stiff equation of state that leads to the formation of an adiabatic core. We impose threshold densities of $10^8$, $10^{10}$, and $10^{12}\,\text{cm}^{-3}$ for the sink particle formation and the onset of the stiff equation of state to study the late, intermediate, and early stages in the evolution of the protostar, respectively. We follow its growth from masses $\simeq 10\,\text{M}_\odot$ to $\simeq 10^5\,\text{M}_\odot$, with an average accretion rate of $\langle\dot{M}_\star\rangle \simeq 2\,\text{M}_\odot\,\text{yr}^{-1}$ for sink particles, and $\simeq 0.8 - 1.4\,\text{M}_\odot\,\text{yr}^{-1}$ for the adiabatic cores. At the end of the simulations, the HII region generated by radiation from the central object has long detached from the protostellar photosphere, but the ionizing radiation remains trapped in the inner host halo, and has thus not yet escaped into the intergalactic medium. Fully coupled, radiation-hydrodynamics simulations hold the key for further progress.

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.

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.

In a previous paper based on dark matter only simulations we show that, in the approximation of an analytic and static potential describing the strongly triaxial and cuspy shape of Milky Way-sized haloes, diffusion due to chaotic mixing in the neighbourhood of the Sun does not efficiently erase phase space signatures of past accretion events. In this second paper we further explore the effect of chaotic mixing using multicomponent Galactic potential models and solar neighbourhood-like volumes extracted from fully cosmological hydrodynamic simulations, thus naturally accounting for the gravitational potential associated with baryonic components, such as the bulge and disc. Despite the strong change in the global Galactic potentials with respect to those obtained in dark matter only simulations, our results confirm that a large fraction of halo particles evolving on chaotic orbits exhibit their chaotic behaviour after periods of time significantly larger than a Hubble time. In addition, significant diffusion in phase space is not observed on those particles that do exhibit chaotic behaviour within a Hubble time.

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 present a model for the evolution of supermassive protostars from their formation at $M_\star \simeq 0.1\,\text{M}_\odot$ until their growth to $M_\star \simeq 10^5\,\text{M}_\odot$. To calculate the initial properties of the object in the optically thick regime we follow two approaches: based on idealized thermodynamic considerations, and on a more detailed one-zone model. Both methods derive a similar value of $n_{\rm F} \simeq 2 \times 10^{17} \,\text{cm}^{-3}$ for the density of the object when opacity becomes important, i.e. the opacity limit. The subsequent evolution of the growing protostar is determined by the accretion of gas onto the object and can be described by a mass-radius relation of the form $R_\star \propto M_\star^{1/3}$ during the early stages, and of the form $R_\star \propto M_\star^{1/2}$ when internal luminosity becomes important. For the case of a supermassive protostar, this implies that the radius of the star grows from $R_\star \simeq 0.65 \,{\rm AU}$ to $R_\star \simeq 250 \,{\rm AU}$ during its evolution. Finally, we use this model to construct a sub-grid recipe for accreting sink particles in numerical simulations. A prime ingredient thereof is a physically motivated prescription for the accretion radius and the effective temperature of the growing protostar embedded inside it. From the latter, we can conclude that photo-ionization feedback can be neglected until very late in the assembly process of the supermassive object.

We present predictions for the quenching of star formation in satellite galaxies of the Local Group from a suite of 30 cosmological zoom simulations of Milky Way-like host galaxies. The Auriga simulations resolve satellites down to the luminosity of the classical dwarf spheroidal galaxies of the Milky Way. We find strong mass-dependent and distance-dependent quenching signals, where dwarf systems beyond 600 kpc are only strongly quenched below a stellar mass of $10^7$ M$_\odot$. Ram pressure stripping appears to be the dominant quenching mechanism and 50% of quenched systems cease star formation within 1 Gyr of first infall. We demonstrate that systems within a host galaxy's $R_{200}$ radius are comprised of two populations: (i) a first infall population that has entered the host halo within the past few Gyrs and (ii) a population of returning `backsplash' systems that have had a much more extended interaction with the host. Backsplash galaxies that do not return to the host galaxy by redshift zero exhibit quenching properties similar to galaxies within $R_{200}$ and are distinct from other external systems. The simulated quenching trend with stellar mass has some tension with observations, but our simulations are able reproduce the range of quenching times measured from resolved stellar populations of Local Group dwarf galaxies.

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 study the incidence of group and filamentary dwarf galaxy accretion into Milky Way (MW) mass haloes using two types of hydrodynamical simulations: EAGLE, which resolves a large cosmological volume, and the AURIGA suite, which are very high resolution zoom-in simulations of individual MW-sized haloes. The present-day 11 most massive satellites are predominantly (75%) accreted in single events, 14% in pairs and 6% in triplets, with higher group multiplicities being unlikely. Group accretion becomes more common for fainter satellites, with 60% of the top 50 satellites accreted singly, 12% in pairs, and 28% in richer groups. A group similar in stellar mass to the Large Magellanic Cloud (LMC) would bring on average 15 members with stellar mass larger than $10^4{~\rm M_\odot}$. Half of the top 11 satellites are accreted along the two richest filaments. The accretion of dwarf galaxies is highly anisotropic, taking place preferentially perpendicular to the halo minor axis, and, within this plane, preferentially along the halo major axis. The satellite entry points tend to be aligned with the present-day central galaxy disc and satellite plane, but to a lesser extent than with the halo shape. Dwarfs accreted in groups or along the richest filament have entry points that show an even larger degree of alignment with the host halo than the full satellite population. We also find that having most satellites accreted as a single group or along a single filament is unlikely to explain the MW disc of satellites.

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.

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.

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.

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.

The stellar disk of the Milky Way shows complex spatial and abundance structure that is central to understanding the key physical mechanisms responsible for shaping our Galaxy. In this study, we use six very high resolution cosmological zoom simulations of Milky Way-sized haloes to study the prevalence and formation of chemically distinct disc components. We find that our simulations develop a clearly bimodal distribution in the $[\rm \alpha/Fe]$ -- $[\rm Fe/H]$ plane. We find two main pathways to creating this dichotomy which operate in different regions of the galaxies: a) an early ($z>1$) and intense high-$\rm[\alpha/Fe]$ star formation phase in the inner region ($R\lesssim 5$ kpc) induced by gas-rich mergers, followed by more quiescent low-$\rm[\alpha/Fe]$ star formation; and b) an early phase of high-$\rm[\alpha/Fe]$ star formation in the outer disc followed by a shrinking of the gas disc owing to a temporarily lowered gas accretion rate, after which disc growth resumes. In process b), a double-peaked star formation history around the time and radius of disc shrinking accentuates the dichotomy. If the early star formation phase is prolonged (rather than short and intense), chemical evolution proceeds as per process a) in the inner region, but the dichotomy is less clear. In the outer region, the dichotomy is only evident if the first intense phase of star formation covers a large enough radial range before disc shrinking occurs; otherwise, the outer disc consists of only low-$\rm[\alpha/Fe]$ sequence stars. We discuss the implication that both processes occurred in the Milky Way.

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 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.

The "Disk of satellites" (DoS) around Milky Way is a highly debated topic with conflicting interpretations of observations and their theoretical models. We perform a comprehensive analysis of all dwarfs detected in the Milky Way and find that the DoS structure depends strongly on the plane identification method and the sample size. In particular, we demonstrate that a smaller sample size produces a higher anisotropy of the spatial distribution and a stronger clustering of the angular momentum of the satellites. Moreover, we calculate the evolution of the 11 classical satellites with proper motion measurements and find that the thin DoS they currently reside in is transient. Furthermore, we analyze two cosmological simulations using the same initial conditions of a Milky Way-sized galaxy, an N-body run with dark matter only and a hydrodynamic one with both baryonic and dark matter, and find that the hydrodynamic simulation produces more anisotropic distribution of satellites than the N-body one. Our results suggest that an anisotropic distribution of satellites in galaxies can originate from baryonic processes in the hierarchical structure formation model, but the claimed highly-flattened, coherently-rotating DoS of the Milky Way may be a small-number selection effect. These findings may help resolve the contradictory claims of DoS in other galaxies and the discrepancy among numerical simulations.

We combine Gaia data release 1 astrometry with Sloan Digital Sky Survey (SDSS) images taken some ~10-15 years earlier, to measure proper motions of stars in the halo of our Galaxy. The SDSS-Gaia proper motions have typical statistical errors of 2 mas/yr down to r ~ 20 mag, and are robust to variations with magnitude and colour. Armed with this exquisite set of halo proper motions, we identify RR Lyrae, blue horizontal branch (BHB), and K giant stars in the halo, and measure their net rotation with respect to the Galactic disc. We find evidence for a gently rotating prograde signal (< $V_\phi$ > ~ 5-25 km/s) in the halo stars, which shows little variation with Galactocentric radius out to 50 kpc. The average rotation signal for the three populations is < $V_\phi$ > = 14 +/- 2 +/- 10 (syst.) km/s. There is also tentative evidence for a kinematic correlation with metallicity, whereby the metal richer BHB and K giant stars have slightly stronger prograde rotation than the metal poorer stars. Using the Auriga simulation suite we find that the old (T >10 Gyr) stars in the simulated halos exhibit mild prograde rotation, with little dependence on radius or metallicity, in general agreement with the observations. The weak halo rotation suggests that the Milky Way has a minor in situ halo component, and has undergone a relatively quiet accretion history.

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.

Observational evidence shows that low-redshift galaxies are surrounded by extended haloes of multiphase gas, the so-called 'circumgalactic medium' (CGM). To study the survival of relatively cool gas (T < 10^5 K) in the CGM, we performed a set of hydrodynamical simulations of cold (T = 10^4 K) neutral gas clouds travelling through a hot (T = 2x10^6 K) and low-density (n = 10^-4 cm^-3) coronal medium, typical of Milky Way-like galaxies at large galactocentric distances (~ 50-150 kpc). We explored the effects of different initial values of relative velocity and radius of the clouds. Our simulations were performed on a two-dimensional grid with constant mesh size (2 pc) and they include radiative cooling, photoionization heating and thermal conduction. We found that for large clouds (radii larger than 250 pc) the cool gas survives for very long time (larger than 250 Myr): despite that they are partially destroyed and fragmented into smaller cloudlets during their trajectory, the total mass of cool gas decreases at very low rates. We found that thermal conduction plays a significant role: its effect is to hinder formation of hydrodynamical instabilities at the cloud-corona interface, keeping the cloud compact and therefore more difficult to destroy. The distribution of column densities extracted from our simulations are compatible with those observed for low-temperature ions (e.g. SiII and SiIII) and for high-temperature ions (OVI) once we take into account that OVI covers much more extended regions than the cool gas and, therefore, it is more likely to be detected along a generic line of sight.

The magnetic fields observed in the Milky~Way and nearby galaxies appear to be in equipartition with the turbulent, thermal, and cosmic ray energy densities, and hence are expected to be dynamically important. However, the origin of these strong magnetic fields is still unclear, and most previous attempts to simulate galaxy formation from cosmological initial conditions have ignored them altogether. Here, we analyse the magnetic fields predicted by the simulations of the Auriga Project, a set of 30 high-resolution cosmological zoom simulations of Milky~Way-like galaxies, carried out with a moving-mesh magneto-hydrodynamics code and a detailed galaxy formation physics model. We find that the magnetic fields grow exponentially at early times owing to a small-scale dynamo with an e-folding time of roughly $100\,\rm{Myr}$ in the center of halos until saturation occurs around $z=2-3$, when the magnetic energy density reaches about $10\%$ of the turbulent energy density with a typical strength of $10-50\,\rm{\mu G}$. In the galactic centers the ratio between magnetic and turbulent energy remains nearly constant until $z=0$. At larger radii, differential rotation in the disks leads to linear amplification that typically saturates around $z=0.5$ to $z=0$. The final radial and vertical variations of the magnetic field strength can be well described by two joint exponential profiles, and are in good agreement with observational constraints. Overall, the magnetic fields have only little effect on the global evolution of the galaxies as it takes too long to reach equipartition. We also demonstrate that our results are well converged with numerical resolution.

We characterize the contribution from accreted material to the galactic discs of the Auriga Project, a set of high resolution magnetohydrodynamic cosmological simulations of late-type galaxies performed with the moving-mesh code AREPO. Our goal is to explore whether a significant accreted (or ex-situ) stellar component in the Milky Way disc could be hidden within the near-circular orbit population, which is strongly dominated by stars born in-situ. One third of our models shows a significant ex-situ disc but this fraction would be larger if constraints on orbital circularity were relaxed. Most of the ex-situ material ($\gtrsim 50\%$) comes from single massive satellites ($> 6 \times 10^{10}~M_{\odot}$). These satellites are accreted with a wide range of infall times and inclination angles (up to $85^{\circ}$). Ex-situ discs are thicker, older and more metal-poor than their in-situ counterparts. They show a flat median age profile, which differs from the negative gradient observed in the in-situ component. As a result, the likelihood of identifying an ex-situ disc in samples of old stars on near-circular orbits increases towards the outskirts of the disc. We show three examples that, in addition to ex-situ discs, have a strongly rotating dark matter component. Interestingly, two of these ex-situ stellar discs show an orbital circularity distribution that is consistent with that of the in-situ disc. Thus, they would not be detected in typical kinematic studies.

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.

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.

It has been suggested that the satellite galaxies of the Milky Way reside in a highly-flattened, kinematically-coherent plane called Disk of Satellites (DoS). The origin of the DoS, however, has been hotly debated, and a number of conflicting claims have been reported in the literature on whether or not the DoS is consistent with predictions from the standard Lambda Cold Dark Matter ($\Lambda$CDM) cosmological model. Here we investigate this issue by comparing a high-resolution, hydrodynamic $\Lambda$CDM simulation of a Milky Way sized galaxy with its dark matter only counterpart. We find the following results: (1) The abundance and distribution of satellite galaxies around a host galaxy is significantly different in the hydro simulation compared to its N-body counterpart; (2) No clear coherent rotation is found in the satellite system, as the fractions of corotating and counter-corotating satellites remain comparable across cosmic time; (3) The satellite distribution evolves significantly with time, from nearly isotropic at high redshift to anisotropic at the present day; (4) The DoS properties strongly depend on sample selection and plane identification methods. Our results imply that the spatially-thin and coherently-rotating DoS reported in Milky Way and other galaxies may be a selection effect of small sample size.

We study the orbital properties of dark matter haloes by combining a spectral method and cosmological simulations of Milky Way-sized galaxies. We compare the dynamics and orbits of individual dark matter particles from both hydrodynamic and $N$-body simulations, and find that the fraction of box, tube and resonant orbits of the dark matter halo decreases significantly due to the effects of baryons. In particular, the central region of the dark matter halo in the hydrodynamic simulation is dominated by regular, short-axis tube orbits, in contrast to the chaotic, box and thin orbits dominant in the $N$-body run. This leads to a more spherical dark matter halo in the hydrodynamic run compared to a prolate one as commonly seen in the $N$-body simulations. Furthermore, by using a kernel based density estimator, we compare the coarse-grained phase-space densities of dark matter haloes in both simulations and find that it is lower by $\sim0.5$ dex in the hydrodynamic run due to changes in the angular momentum distribution, which indicates that the baryonic process that affects the dark matter is irreversible. Our results imply that baryons play an important role in determining the shape, kinematics and phase-space density of dark matter haloes in galaxies.

We analyse the properties of the HI gas distribution in the Auriga project, a set of magnetohydrodynamic cosmological simulations performed with the moving-mesh code AREPO and a physics model for galaxy formation that succeeds in forming realistic late-type galaxies in the 30 Milky Way-sized haloes simulated in this project. We use a simple approach to estimate the neutral hydrogen fraction in our simulation set, which treats low-density and star-forming gas separately, and we explore two different prescriptions to subtract the contribution of molecular hydrogen from the total HI content. The HI gas in the vast majority of the systems forms extended discs although more disturbed morphologies are present. Notwithstanding the general good agreement with observed HI properties -- such as radial profiles and the mass-diameter relation -- the Auriga galaxies are systematically larger and more gas-rich than typical nearby galaxies. Interestingly, the amount of HI gas outside the disc plane correlates with the star formation rate, consistent with a picture where most of this extra-planar HI gas originates from a fountain-like flow. Our findings are robust with respect to the different assumptions adopted for computing the molecular hydrogen fraction and do not vary significantly over a wide range of numerical resolution. The HI modelling introduced in this paper can be used in future work to build artificial interferometric HI data cubes, allowing an even closer comparison of the gas dynamics in simulated galaxies with observations.

We study the impact of baryons on the distribution of dark matter in a Milky Way-size halo by comparing a high-resolution, moving-mesh cosmological simulation with its dark matter-only counterpart. We identify three main processes related to baryons -- adiabatic contraction, tidal disruption and reionization -- which jointly shape the dark matter distribution in both the main halo and its subhalos. The relative effect of each baryonic process depends strongly on the subhalo mass. For massive subhalos with maximum circular velocity $v_{\rm max} > 35 km/s$, adiabatic contraction increases the dark matter concentration, making these halos less susceptible to tidal disruption. For low-mass subhalos with $v_{\rm max} < 20 km/s$, reionization effectively reduces their mass on average by $\approx$ 30% and $v_{\rm max}$ by $\approx$ 20%. For intermediate subhalos with $20 km/s < v_{\rm max} < 35 km/s$, which share a similar mass range as the classical dwarf spheroidals, strong tidal truncation induced by the main galaxy reduces their $v_{\rm max}$. Moreover, the stellar disk of the main galaxy effectively depletes subhalos near the central region. As a combined result of reionization and increased tidal disruption, the total number of low-mass subhalos in the hydrodynamic simulation is nearly halved compared to that of the $\textit{N-}$body simulation. We do not find dark matter cores in dwarf galaxies, unlike previous studies that employed bursty feedback-driven outflows. The substantial impact of baryons on the abundance and internal structure of subhalos suggests that galaxy formation and evolution models based on $\textit{N}$-body simulations should include these physical processes as major components.

We introduce a suite of thirty cosmological magneto-hydrodynamical zoom simulations of the formation of galaxies in isolated Milky Way mass dark haloes. These were carried out with the moving mesh code \textlcsc{AREPO}, together with a comprehensive model for galaxy formation physics, including AGN feedback and magnetic fields, which produces realistic galaxy populations in large cosmological simulations. We demonstrate that our simulations reproduce a wide range of present-day observables, in particular, two component disc dominated galaxies with appropriate stellar masses, sizes, rotation curves, star formation rates and metallicities. We investigate the driving mechanisms that set present-day disc sizes/scale lengths, and find that they are related to the angular momentum of halo material. We show that the largest discs are produced by quiescent mergers that inspiral into the galaxy and deposit high angular momentum material into the pre-existing disc, simultaneously increasing the spin of dark matter and gas in the halo. More violent mergers and strong AGN feedback play roles in limiting disc size by destroying pre-existing discs and by suppressing gas accretion onto the outer disc, respectively. The most important factor that leads to compact discs, however, is simply a low angular momentum for the halo. In these cases, AGN feedback plays an important role in limiting central star formation and the formation of a massive bulge.

Recent studies have revealed an oscillating asymmetry in the vertical structure of the Milky Way's disc. Here we analyze 16 high-resolution, fully cosmological simulations of the evolution of individual Milky Way-sized galaxies, carried out with the MHD code AREPO. At redshift zero, about $70\%$ of our galactic discs show strong vertical patterns, with amplitudes that can exceed 2 kpc. Half of these are typical `integral sign' warps. The rest are oscillations similar to those observed in the Milky Way. Such structures are thus expected to be common. The associated mean vertical motions can be as large as 30~km/s. Cold disc gas typically follows the vertical patterns seen in the stars. These perturbations have a variety of causes: close encounters with satellites, distant flybys of massive objects, accretion of misaligned cold gas from halo infall or from mergers. Tidally induced vertical patterns can be identified in both young and old stellar populations, whereas those originating from cold gas accretion are seen mainly in the younger populations. Galaxies with regular or at most weakly perturbed discs are usually, but not always, free from recent interactions with massive companions, although we have one case where an equilibrium compact disc reforms after a merger.

Numerical simulations of structure formation have recorded a remarkable progress in the recent years, in particular due to the inclusion of baryonic physics evolving with the dark matter component. We generate Monte Carlo realizations of the dark matter sub-halo population based on the results of the recent hydrodynamical simulation suite of Milky Way-sized galaxies. We then simulate the gamma-ray sky for both the setup of the 3FGL and 2FHL Fermi Large Area Telescope (LAT) catalogs, including the contribution from the annihilation of dark matter in the sub-halos. We find that the flux sensitivity threshold strongly depends on the particle dark matter mass, and more mildly also on its annihilation channel and the observation latitude. The results differ for the 3FGL and 2FHL catalogs, given their different energy thresholds. We also predict that the number of dark matter sub-halos among the unassociated sources is very small. A null number of detectable sub-halos in the Fermi-LAT 3FGL catalog would imply upper limits on the dark matter annihilation cross section into $b\bar{b}$ of $2 \cdot 10^{-26}$ ($5 \cdot 10^{-25}$) cm$^3$/s with $M_{\rm DM}$= 50 (1000) GeV. We find less than one extended sub-halo in the Fermi-LAT 3FGL catalog. As a matter of fact, the differences in the spatial and mass distribution of sub-halos between hydrodynamic and dark matter-only runs do not have significant impact on the gamma-ray dark matter phenomenology.

In star-forming galaxies, stellar feedback can have a dual effect on the circumgalactic medium both suppressing and stimulating gas accretion. The trigger of gas accretion can be caused by disc material ejected into the halo in the form of fountain clouds and by its interaction with the surrounding hot corona. Indeed, at the disc-corona interface, the mixing between the cold/metal-rich disc gas (T <~ 10^4 K) and the hot coronal gas (T >~ 10^6 K) can dramatically reduce the cooling time of a portion of the corona and produce its condensation and accretion. We studied the interaction between fountain clouds and corona in different galactic environments through parsec-scale hydrodynamical simulations, including the presence of thermal conduction, a key mechanism that influences gas condensation. Our simulations showed that the coronal gas condensation strongly depends on the galactic environment, in particular it is less efficient for increasing virial temperature/mass of the haloes where galaxies reside and it is fully ineffective for objects with virial masses larger than 10^13 Msun. This result implies that the coronal gas cools down quickly in haloes with low-intermediate virial mass (Mvir <~ 3 x 10^12 Msun) but the ability to cool the corona decreases going from late-type to early-type disc galaxies, potentially leading to the switching off of accretion and the quenching of star formation in massive systems.

We study the formation of galactic outflows from supernova explosions (SNe) with the moving-mesh code AREPO in a stratified column of gas with a surface density similar to the Milky Way disk at the solar circle. We compare different simulation models for SNe placement and energy feedback, including cosmic rays (CR), and find that models that place SNe in dense gas and account for CR diffusion are able to drive outflows with similar mass loading as obtained from a random placement of SNe with no CRs. Despite this similarity, CR-driven outflows differ in several other key properties including their overall clumpiness and velocity. Moreover, the forces driving these outflows originate in different sources of pressure, with the CR diffusion model relying on non-thermal pressure gradients to create an outflow driven by internal pressure and the random-placement model depending on kinetic pressure gradients to propel a ballistic outflow. CRs therefore appear to be non-negligible physics in the formation of outflows from the interstellar medium.

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.

We study the vertical structure of a stellar disk obtained from a fully cosmological high-resolution hydrodynamical simulation of the formation of a Milky Way-like galaxy. At the present day, the disk's mean vertical height shows a well-defined and strong pattern, with amplitudes as large as 3 kpc in its outer regions. This pattern is the result of a satellite - host halo - disk interaction and reproduces, qualitatively, many of the observable properties of the Monoceros Ring. In particular we find disk material at the distance of Monoceros ($R \sim$ 12-16 kpc, galactocentric) extending far above the mid plane ( 30$^{\circ}$, $\langle Z \rangle \sim$ 1-2 kpc) in both hemispheres, as well as well-defined arcs of disk material at heliocentric distances $\gtrsim 5$ kpc. The pattern was first excited $\approx 3$ Gyr ago as an $m=1$ mode that later winds up into a leading spiral pattern. Interestingly, the main driver behind this perturbation is a low-mass low-velocity fly-by encounter. The satellite has total mass, pericentre distance and pericentric velocity of $\sim 5\%$ of the host, $\sim 80$ kpc, and 215 km/s, respectively. The satellite is not massive enough to directly perturb the galactic disk but we show that the density field of the host dark matter halo responds to this interaction resulting in a strong amplification of the perturbative effects. This subsequently causes the onset and development of the Monoceros-like feature.

We use a high resolution cosmological zoom simulation of a Milky Way-sized halo to study the observable features in velocity and metallicity space associated with the dynamical influence of spiral arms. For the first time, we demonstrate that spiral arms, that form in a disc in a fully cosmological environment with realistic galaxy formation physics, drive large-scale systematic streaming motions. In particular, on the trailing edge of the spiral arms the peculiar galacto-centric radial and azimuthal velocity field is directed radially outward and azimuthally backward, whereas it is radially inward and azimuthally forward on the leading edge. Owing to the negative radial metallicity gradient, this systematic motion drives, at a given radius, an azimuthal variation in the residual metallicity that is characterised by a metal rich trailing edge and a metal poor leading edge. We show that these signatures are theoretically observable in external galaxies with Integral Field Unit instruments such as VLT/MUSE, and if detected, would provide evidence for large-scale systematic radial migration driven by spiral arms.

Vertically extended, high velocity dispersion stellar distributions appear to be a ubiquitous feature of disc galaxies, and both internal and external mechanisms have been proposed to be the major driver of their formation. However, it is unclear to what extent each mechanism can generate such a distribution, which is likely to depend on the assembly history of the galaxy. To this end, we perform 16 high resolution cosmological-zoom simulations of Milky Way-sized galaxies using the state-of-the-art cosmological magneto-hydrodynamical code \textlcsc{AREPO}, and analyse the evolution of the vertical kinematics of the stellar disc in connection with various heating mechanisms. We find that the bar is the dominant heating mechanism in most cases, whereas spiral arms, radial migration, and adiabatic heating from mid-plane density growth are all sub-dominant. The strongest source, though less prevalent than bars, originates from external perturbations from satellites/sub-halos of masses log$_{10} (M/\rm M_{\odot}) \gtrsim 10$. However, in many simulations the orbits of newborn star particles become cooler with time, such that they dominate the shape of the age-velocity dispersion relation and overall vertical disc structure unless a strong external perturbation takes place.

A recent observational study of haloes of nearby Milky Way-like galaxies shows that only half (four out of eight) of the current sample exhibits strong negative metallicity ([Fe/H]) gradients. This is at odds with predictions from hydrodynamical simulations where such gradients are ubiquitous. In this Letter, we use high resolution cosmological hydrodynamical simulations to study the [Fe/H] distribution of galactic haloes. We find that kinematically selected stellar haloes, including both in-situ and accreted particles, have an oblate [Fe/H] distribution. Spherical [Fe/H] radial profiles show strong negative gradients within 100 kpc, in agreement with previous numerical results. However, the projected median [Fe/H] profiles along the galactic disc minor axis, typically obtained in observations, are significantly flatter. The median [Fe/H] values at a given radius are larger for the spherical profiles than for the minor axis profiles by as much as 0.4 dex within the inner 50 kpc. Similar results are obtained if only the accreted stellar component is considered indicating that the differences between spherical and minor axis profiles are not purely driven by heated disc star particles formed in situ. Our study highlights the importance of performing careful comparisons between models and observations of halo [Fe/H] distributions.

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.

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.

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.

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 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.

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.

High-velocity clouds consist of cold gas that appears to be raining down from the halo to the disc of the Milky Way. Over the past fifty years, two competing scenarios have attributed their origin either to gas accretion from outside the Galaxy or to circulation of gas from the Galactic disc powered by supernova feedback (galactic fountain). Here we show that both mechanisms are simultaneously at work. We use a new galactic fountain model combined with high-resolution hydrodynamical simulations. We focus on the prototypical cloud complex C and show that it was produced by an explosion that occurred in the Cygnus-Outer spiral arm about 150 million years ago. The ejected material has triggered the condensation of a large portion of the circumgalactic medium and caused its subsequent accretion onto the disc. This fountain-driven cooling of the lower Galactic corona provides the low-metallicity gas required by chemical evolution models of the Milky Way's disc.

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 analyse the properties of the circum-galactic medium and the metal content of the stars comprising the central galaxy in eight hydrodynamical `zoom-in' simulations of disc galaxy formation. We use these properties as a benchmark for our model of galaxy formation physics implemented in the moving-mesh code AREPO, which succeeds in forming quite realistic late-type spirals in the set of `Aquarius' initial conditions of Milky Way-sized haloes. Galactic winds significantly influence the morphology of the circum-galactic medium and induce bipolar features in the distribution of heavy elements. They also affect the thermodynamic properties of the circum-galactic gas by supplying an energy input that sustains its radiative losses. Although a significant fraction of the heavy elements are transferred from the central galaxy to the halo, and even beyond the virial radius, enough metals are retained by stars to yield a peak in their metallicity distributions at about $Z_{\odot}$. All our default runs overestimate the stellar [O/Fe] ratio, an effect that we demonstrate can be rectified by an increase of the adopted SN type Ia rate. Nevertheless, the models have difficulty in producing stellar metallicity gradients of the same strength as observed in the Milky Way.

The relative impact of radiation pressure and photoionization feedback from young stars on surrounding gas is studied with hydrodynamic radiative transfer (RT) simulations. The calculations focus on the single-scattering (direct radiation pressure) and optically thick regime, and adopt a moment-based RT-method implemented in the moving-mesh code AREPO. The source luminosity, gas density profile and initial temperature are varied. At typical temperatures and densities of molecular clouds, radiation pressure drives velocities of order ~20 km/s over 1-5 Myr; enough to unbind the smaller clouds. However, these estimates ignore the effects of photoionization that naturally occur concurrently. When radiation pressure and photoionization act together, the latter is substantially more efficient, inducing velocities comparable to the sound speed of the hot ionized medium (10-15 km/s) on timescales far shorter than required for accumulating similar momentum with radiation pressure. This mismatch allows photoionization to dominate the feedback as the heating and expansion of gas lowers the central densities, further diminishing the impact of radiation pressure. Our results indicate that a proper treatment of the impact of young stars on the interstellar medium needs to primarily account for their ionization power whereas direct radiation pressure appears to be a secondary effect. This conclusion may change if extreme boosts of the radiation pressure by photon trapping are assumed.

Observationally, magnetic fields reach equipartition with thermal energy and cosmic rays in the interstellar medium of disk galaxies such as the Milky Way. However, thus far cosmological simulations of the formation and evolution of galaxies have usually neglected magnetic fields. We employ the moving-mesh code \textsc{Arepo} to follow for the first time the formation and evolution of a Milky Way-like disk galaxy in its full cosmological context while taking into account magnetic fields. We find that a prescribed tiny magnetic seed field grows exponentially by a small-scale dynamo until it saturates around $z=4$ with a magnetic energy of about $10\%$ of the kinetic energy in the center of the galaxy's main progenitor halo. By $z=2$, a well-defined gaseous disk forms in which the magnetic field is further amplified by differential rotation, until it saturates at an average field strength of $\sim 6 \mug$ in the disk plane. In this phase, the magnetic field is transformed from a chaotic small-scale field to an ordered large-scale field coherent on scales comparable to the disk radius. The final magnetic field strength, its radial profile and the stellar structure of the disk compare well with observational data. A minor merger temporarily increases the magnetic field strength by about a factor of two, before it quickly decays back to its saturation value. Our results are highly insensitive to the initial seed field strength and suggest that the large-scale magnetic field in spiral galaxies can be explained as a result of the cosmic structure formation process.

We present cosmological hydrodynamical simulations of eight Milky Way-sized haloes that have been previously studied with dark matter only in the Aquarius project. For the first time, we employ the moving-mesh code AREPO in zoom simulations combined with a comprehensive model for galaxy formation physics designed for large0 cosmological simulations. Our simulations form in most of the eight haloes strongly disc-dominated systems with realistic rotation curves, close to exponential surface density profiles, a stellar-mass to halo-mass ratio that matches expectations from abundance matching techniques, and galaxy sizes and ages consistent with expectations from large galaxy surveys in the local Universe. There is no evidence for any dark matter core formation in our simulations, even so they include repeated baryonic outflows by supernova-driven winds and black hole quasar feedback. For one of our haloes, the object studied in the recent `Aquila' code comparison project, we carried out a resolution study with our techniques, covering a dynamic range of 64 in mass resolution. Without any change in our feedback parameters, the final galaxy properties are reassuringly similar, in contrast to other modelling techniques used in the field that are inherently resolution dependent. This success in producing realistic disc galaxies is reached, in the context of our interstellar medium treatment, without resorting to a high density threshold for star formation, a low star formation efficiency, or early stellar feedback, factors deemed crucial for disc formation by other recent numerical studies.

The spatial segregation between dSphs and dIrrs in the Local Group has long been regarded as evidence of an interaction with their host galaxies. In this paper, we assume that ram-pressure stripping is the dominant mechanism that removed gas from the dSphs and we use this to derive a lower bound on the density of the corona of the Milky Way at large distances (50-90 kpc) from the Galactic centre. At the same time, we derive an upper bound by demanding that the interstellar medium of the dSphs is in pressure equilibrium with the hot corona. We consider two dwarfs (Sextans and Carina) with well-determined orbits and star formation histories. Our approach introduces several novel features: we use the measured star formation histories of the dwarfs to derive the time at which they last lost their gas, and (via a modified version of the Kennicutt-Schmidt relation) their internal gas density at that time; we use a large suite of 2D hydrodynamical simulations to model the gas stripping; and we include supernova feedback tied to the gas content. Despite having very different orbits and star formation histories, we find results for the two dSphs that are in excellent agreement with one another. We derive an average particle density of the corona of the Milky Way at 50-90 kpc in the range 1.3-3.6 10^{-4} cm^{-3}. Including additional constraints from X-ray emission limits and pulsar dispersion measurements, we extrapolate Galactic coronal density profiles and we estimate the fraction of baryons that can exist within the virial radius of the Milky Way. For an isothermal corona (T=1.8 10^6 K) this is small, 10-20 % of the universal baryon fraction. Only a hot (T=3 10^6 K) and adiabatic corona can contain all of the Galaxy's missing baryons. Models for the Milky Way must explain why its corona is in a hot adiabatic thermal state or why a large fraction of its baryons lie beyond the virial radius.

Disc galaxies like the Milky Way are expected to be surrounded by massive coronae of hot plasma that may contain a significant fraction of the so-called missing baryons. We investigate whether the local (|vLSR|<400 km/s) warm-hot absorption features observed towards extra-Galactic sources or halo stars are consistent with being produced by the cooling of the Milky Way's corona. In our scheme, cooling occurs at the interface between the disc and the corona and it is triggered by positive supernova feedback. We combine hydrodynamical simulations with a dynamical 3D model of the galactic fountain to predict the all-sky distribution of this cooling material, and we compare it with the observed distribution of detections for different `warm' (SiIII, SiIV, CII, CIV) and `hot' (OVI) ionised species. The model reproduces the position-velocity distribution and the column densities of the vast majority of warm absorbers and about half of OVI absorbers. We conclude that the warm-hot gas responsible for most of the detections lies within a few kiloparsecs from the Galactic plane, where high-metallicity material from the disc mixes efficiently with the hot corona. This process provides an accretion of a few Mo/yr of fresh gas that can easily feed the star formation in the disc of the Galaxy. The remaining OVI detections are likely to be a different population of absorbers, located in the outskirts of the Galactic corona and/or in the circumgalactic medium of nearby galaxies.

We show that the ultraviolet absorption features, newly discovered in HST spectra, are consistent with being formed in a layer that extends a few kpc above the disk of the Milky Way. In this interface between the disk and the Galactic corona, high-metallicity gas ejected from the disk by supernova feedback can mix efficiently with the virial-temperature coronal material. The mixing process triggers the cooling of the lower corona down to temperatures encompassing the characteristic range of the observed absorption features, producing a net supernova-driven gas accretion onto the disk at a rate of a few Msun/yr. We speculate that this mechanism explains how the hot-mode of cosmological accretion feeds star formation in galactic disks.

Accretion of fresh gas at a rate of ~ 1 M_{sun} yr^{-1} is necessary in star-forming disc galaxies, such as the Milky Way, in order to sustain their star-formation rates. In this work we present the results of a new hydrodynamic simulation supporting the scenario in which the gas required for star formation is drawn from the hot corona that surrounds the star-forming disc. In particular, the cooling of this hot gas and its accretion on to the disc are caused by the passage of cold galactic fountain clouds through the corona.

In galaxies like the Milky Way, cold (~ 10^4 K) gas ejected from the disc by stellar activity (the so-called galactic-fountain gas) is expected to interact with the virial-temperature (~ 10^6 K) gas of the corona. The associated transfer of momentum between cold and hot gas has important consequences for the dynamics of both gas phases. We quantify the effects of such an interaction using hydrodynamical simulations of cold clouds travelling through a hot medium at different relative velocities. Our main finding is that there is a velocity threshold between clouds and corona, of about 75 km/s, below which the hot gas ceases to absorb momentum from the cold clouds. It follows that in a disc galaxy like the Milky Way a static corona would be rapidly accelerated: the corona is expected to rotate and to lag, in the inner regions, by ~ 80-120 km/s with respect to the cold disc. We also show how the existence of this velocity threshold can explain the observed kinematics of the cold extra-planar gas.

It is argued that galaxies like ours sustain their star formation by transferring gas from an extensive corona to the star-forming disc. The transfer is effected by the galactic fountain -- cool clouds that are shot up from the plane to kiloparsec heights above the plane. The Kelvin-Helmholtz instability strips gas from these clouds. If the pressure and the the metallicity of the corona are high enough, the stripped gas causes a similar mass of coronal gas to condense in the cloud's wake. Hydrodynamical simulations of cloud-corona interaction are presented. These confirm the existence of a critical ablation rate above which the corona is condensed, and imply that for the likely parameters of the Galactic corona this rate lies near the actual ablation rate of clouds. In external galaxies trails of HI behind individual clouds will not be detectable, although the integrated emission from all such trails should be significant. Parts of the trails of the clouds that make up the Galaxy's fountain should be observable and may account for features in targeted 21-cm observations of individual high-velocity clouds and surveys of Galactic HI emission. Taken in conjunction with the known decline in the availability of cold infall with increasing cosmic time and halo mass, the proposed mechanism offers a promising explanation of the division of galaxies between the blue cloud to the red sequence in the colour-luminosity plane.

Star-forming disc galaxies such as the Milky Way need to accrete $\gsim$ 1 $M_{\odot}$ of gas each year to sustain their star formation. This gas accretion is likely to come from the cooling of the hot corona, however it is still not clear how this process can take place. We present simulations supporting the idea that this cooling and the subsequent accretion are caused by the passage of cold galactic-fountain clouds through the hot corona. The Kelvin-Helmholtz instability strips gas from these clouds and the stripped gas causes coronal gas to condense in the cloud's wake. For likely parameters of the Galactic corona and of typical fountain clouds we obtain a global accretion rate of the order of that required to feed the star formation.

The kinematics of the extra-planar neutral and ionised gas in disc galaxies shows a systematic decline of the rotational velocity with height from the plane (vertical gradient). This feature is not expected for a barotropic gas, whilst it is well reproduced by baroclinic fluid homogeneous models. The problem with the latter is that they require gas temperatures (above $10^5$ K) much higher than the temperatures of the cold and warm components of the extra-planar gas layer. In this paper, we attempt to overcome this problem by describing the extra-planar gas as a system of gas clouds obeying the Jeans equations. In particular, we consider models having the observed extra-planar gas distribution and gravitational potential of the disc galaxy NGC 891: for each model we construct pseudo-data cubes and we compare them with the HI data cube of NGC 891. In all cases the rotational velocity gradients are in qualitative agreement with the observations, but the synthetic and the observed data cubes of NGC 891 show systematic differences that cannot be accommodated by any of the explored models. We conclude that the extra-planar gas in disc galaxies cannot be satisfactorily described by a stationary Jeans-like system of gas clouds.

In a previous paper the complex-shift method has been applied to self-gravitating spherical systems, producing new analytical axisymmetric density-potential pairs. We now extend the treatement to the Miyamoto-Nagai disc and to the Binney logarithmic halo, and we study the resulting axisymmetric and triaxial analytical density-potential pairs; we also show how to obtain the surface density of shifted systems from the complex-shift of the surface density of the parent model. In particular, the systems obtained from Miyamoto-Nagai discs can be used to describe disc galaxies with a peanut-shaped bulge or with a central triaxial bar, depending on the direction of the shift vector. By using a constructive method that can be applied to generic axisymmetric systems, we finally show that the Miyamoto-Nagai and the Satoh discs, and the Binney logarithmic halo, cannot be obtained from the complex-shift of any spherical parent distribution. As a by-product of this study we also found two new generating functions in closed form for even and odd Legendre polynomials, respectively.