The MillenniumTNG Project: the hydrodynamical full physics simulation and a first look at its galaxy clusters

ABSTRACT

Cosmological simulations are an important theoretical pillar for understanding non-linear structure formation in our Universe and for relating it to observations on large scales. In several papers, we introduce our MillenniumTNG (MTNG) project that provides a comprehensive set of high-resolution, large-volume simulations of cosmic structure formation aiming to better understand physical processes on large scales and to help interpret upcoming large-scale galaxy surveys. We here focus on the full physics box MTNG740 that computes a volume of |$740\, \mathrm{Mpc}^3$| with a baryonic mass resolution of |$3.1\times ~10^7\, \mathrm{M_\odot }$| using arepo with 80.6 billion cells and the IllustrisTNG galaxy formation model. We verify that the galaxy properties produced by MTNG740 are consistent with the TNG simulations, including more recent observations. We focus on galaxy clusters and analyse cluster scaling relations and radial profiles. We show that both are broadly consistent with various observational constraints. We demonstrate that the SZ-signal on a deep light-cone is consistent with Planck limits. Finally, we compare MTNG740 clusters with galaxy clusters found in Planck and the SDSS-8 RedMaPPer richness catalogue in observational space, finding very good agreement as well. However, simultaneously matching cluster masses, richness, and Compton-y requires us to assume that the SZ mass estimates for Planck clusters are underestimated by 0.2 dex on average. Due to its unprecedented volume for a high-resolution hydrodynamical calculation, the MTNG740 simulation offers rich possibilities to study baryons in galaxies, galaxy clusters, and in large-scale structure, and in particular their impact on upcoming large cosmological surveys.

1 INTRODUCTION

Recent large cosmological hydrodynamical simulations have been very successful in reproducing realistic galaxy populations on cosmological scales. These include Illustris (Vogelsberger et al. 2014a, b), Eagle (Schaye et al. 2015), HorizonAGN (Dubois et al. 2016), IllustrisTNG (Springel et al. 2018), Simba (Davé et al. 2019), NewHorizon (Dubois et al. 2021), Thesan (Kannan et al. 2022b), and Astrid (Bird et al. 2022; Ni et al. 2022). (For a recent overview of the corresponding modelling techniques, see Vogelsberger et al. 2020.) However, among these simulation boxes even TNG300 – as the largest of them – is still too small to properly study large-scale structure features such as baryonic acoustic oscillations or to contain a representative sample of massive galaxy clusters. Larger cosmological simulations exist, e.g. Magneticum (Dolag, Komatsu & Sunyaev 2016) or Bahamas (McCarthy et al. 2017), but they can only afford a mass resolution that is too low to properly resolve individual galaxies, or have not reached z = 0 yet (Astrid). Yet, it is crucial to understand baryonic physics on large scales, i.e. in galaxy clusters or in the cosmic web, to anchor the baryonic physics model with galaxies and at low redshift where we have many good observational constraints.

Moreover, even larger cosmological box simulations are needed as a basis for the interpretation of upcoming enormous galaxy surveys which aim to better constrain the cosmological parameters of the Universe to the per cent level (e.g. DES, eBOSS, DESI, or Euclid). However, such simulations with sufficiently large volumes can only follow dark matter (Potter, Stadel & Teyssier 2017; Angulo & Hahn 2022) and require significant post-processing to model actual observables, which for the most part are either exclusively based on baryons or are at least affected by baryonic physics.

A similar problem arises when simulating galaxy clusters. Their sparsity requires large simulation volumes in order to contain a significant number of massive clusters. To simulate them at sufficient resolution, zoom simulations are thus often employed that focus most resolution elements in a small region of interest centred around a single galaxy cluster. These zoom simulations (Bahé et al. 2017; Barnes et al. 2017; Cui et al. 2018) enable the study of the internal structure of galaxy clusters and galaxy cluster scaling relations without the need to follow large volumes at high resolution. It is difficult, however, to compose truly representative samples of clusters with this technique, and to properly model their cosmological foregrounds at the same time.

This study is one of the introductory papers of our new MillenniumTNG project in which we seek to make progress on pushing direct hydrodynamical simulations of galaxy formation to much larger volume than available thus far, and on linking these hydrodynamical results to still larger dark-matter-only simulations, thereby allowing the hydrodynamical results to be more reliably extrapolated to cosmological scales. Simultaneously, this offers a new opportunity to study large samples of hydrodynamically simulated galaxy clusters.

To this end, our work applies the IllustrisTNG state-of-the-art cosmological galaxy formation model in an unparalleled large simulation volume at a mass resolution that still allows us to reasonably describe the properties of individual galaxies. Our most ambitious hydrodynamic calculation is carried out in the |$500\, h^{-1}\, \mathrm{Mpc}$| periodic box size of the seminal Millennium simulation, thus motivating the name ‘MillenniumTNG’ we coined for the whole project. The nearly 15 times increase of the simulated volume compared to TNG300 at slightly lower mass resolution allows our hydrodynamic calculation to be more directly compared to upcoming large volume cosmological surveys, and importantly, it enables us to calibrate and improve approximate methods to predict galaxy catalogues based on dark-matter-only simulations. To facilitate the latter, MillenniumTNG additionally consists of a suite of dark-matter-only simulations computed in the same volume, with two ‘fixed-and-paired’ versions at each resolution that use a variance suppression technique that effectively boosts the statistical power of the simulated volume (Angulo & Pontzen 2016). Furthermore, we have computed yet much larger dark-matter-only models that also include massive neutrinos, with up to 1.1 trillion particles in a |$3000\, \mathrm{Mpc}$| box, which is meant to propel the statistical power of our predictions into the regime probed by the upcoming surveys.

We refer to a companion paper by Hernández-Aguayo et al. (2022) for full details on the simulation suite, its data products (including light-cone outputs and merger trees) and a study of basic matter and halo clustering statistics. In this paper, we focus on introducing the flagship hydrodynamical full physics simulation of the MillenniumTNG project, and on giving an initial characterization of the galaxy clusters in the simulation. This offers a first glimpse at the possibilities the project offers to understand current and future cosmological observations. In Barrera et al. (2022), we present a novel version of the L-Galaxies semi-analytic model of galaxy formation and its application to the light-cone outputs of the MTNG simulations. In further companion papers, Kannan et al. (2022a) analyse the properties of very high-redshift galaxies, while Bose et al. (2022) present a galaxy clustering study based on colour-selected galaxy samples. Hadzhiyska et al. (2022a, b) examine aspects of galaxy assembly bias, whereas Delgado et al. (in preparation) study intrinsic alignments and galaxy shapes. Contreras et al. (2022) introduce an inference technique to constrain the cosmological parameters of MTNG from galaxy clustering. Finally, Ferlito et al. (in preparation) study weak gravitational lensing both in the dark matter and the full physics runs.

This paper is structured as follows. Section 2 introduces the MillenniumTNG full physics box and summarizes its galaxy formation model and model parameters. Section 3 compares global galaxy properties of the MTNG740 full physics box to the TNG simulations, as well as to updated observational data. Section 4 discusses galaxy clusters in MTNG740 and presents various galaxy cluster scaling relations and compares them to observations. Section 5 analyses different radial profiles of galaxy clusters and puts them in relation to observations. Section 6 compares the SZ signal of a galaxy cluster at z = 0.25 computed from a snapshot time-slice with spherical and cylindrical apertures to the same observable but calculated using the simulation data on a deep full backwards light-cone. Section 7 considers galaxy and cluster observables to compare the galaxy clusters in MTNG740 in observational space with clusters observed with Planck and SDSS-8. Finally, Section 8 concludes with a summary of the paper and an outlook on applications of the MTNG simulations.

2 FULL PHYSICS SIMULATIONS IN MILLENNIUM-TNG

The MTNG physics model is based on the IllustrisTNG (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019b; Pillepich et al. 2018b, 2019; Springel et al. 2018) galaxy formation model (Weinberger et al. 2017; Pillepich et al. 2018a) that has been shown to produce a fairly realistic galaxy population on cosmological scales. It includes primordial and metal line cooling (Vogelsberger et al. 2013), an explicit sub-grid model for the interstellar medium and star formation (Springel & Hernquist 2003), mass return from stars and metal enrichment of the insterstellar medium by core-collapse supernovae, thermonuclear supernovae, and AGB stars, an effective model for galactic winds (Pillepich et al. 2018a), and a model for the creation and growth of supermassive black holes as well as their feedback as active galactic nuclei (AGNs; Weinberger et al. 2017).

We keep the physics model of IllustrisTNG with all its parameters unchanged, subject to only a few minor modifications. We fixed all the small issues found during and after the IllustrisTNG simulations have been run and that are documented in Nelson et al. (2019a). Notably, we corrected the unintentional abrupt start of the UV background radiation at the epoch of cosmic reionization. We needed to remove magnetic fields and individual metal species from the model to reduce the memory requirements of the simulation. Instead of following individual metal species we reverted to evolving only a single scalar field that tracks the total metallicity of cells and star particles. Based on small test runs, we do not expect this change to affect the simulation results for galaxies in any significant way. Nevertheless, removing the magnetic fields does however change the model slightly (Pillepich et al. 2018a), and we discuss related effects in more detail in Section 3. Since magnetic fields are omitted, we switched to using an exact Riemann solver for hydrodynamics instead of using an approximate one. In contrast to IllustrisTNG, we did not include any passive tracer particles, again for reasons of memory consumption.

Despite keeping the physics model almost identical to IllustrisTNG, we implemented a large number of technical changes and improvements to the arepo code (Springel 2010; Pakmor et al. 2016; Weinberger, Springel & Pakmor 2020) in order to make the MTNG740 full physics box fit into the memory of the supercomputer available to us. For example, for this purpose, we now use shared memory on compute nodes via MPI-3 to store identical data only once per node rather than once per compute core. This includes, most importantly, information about the domain decomposition and the top-level tree that is shared between all MPI ranks, as well as various data tables such as stellar yields or photometric tables for stellar population synthesis models. We now also exploit the shared memory on nodes to improve the efficiency of various global MPI operations, for example, by computing collective results first on each node, and then exchanging them only between nodes rather than between all cores individually. Similarly, we have replaced the global domain decomposition of the code with a hierarchical strategy that first subdivides the simulated volume among all nodes, and then cuts down the pieces further within each node to distribute the work among the available compute cores. Despite the use of well over 10⁵ MPI ranks, this allows us to always efficiently find a domain decomposition with a well-balanced computational load and a maximum memory overhead of around |$10{{\ \rm per\ cent}}$|⁠. Finally, we have reordered some operations in the gravitational tree algorithm, allowing the front-loading of a communication step. In this way, we can now almost completely avoid wait times originating from different numbers of imported and exported resolution elements on a local MPI task (Springel2010).

We have included as much post-processing as possible already on the fly during the simulation. This includes running the friends-of-friends (FoF) group finder and a novel SUBFIND-HBT substructure finder (Springel et al. 2021) adopted from the Gadget4 into our moving-mesh code Arepo. It also involves computing merger trees and matter power spectra already while the simulation is run. Finally, we output resolution elements when they intersect the past backwards light-cone of a fiducial observer, i.e. light emitted from them will be seen by the fiducial observer exactly at z = 0. The box is periodically replicated if the light-cone extends beyond the size of the box. We use five different light-cone geometries as described in detail in Hernández-Aguayo et al. (2022).

In the following text, we focus almost exclusively on the hydrodynamical full physics box of the MillenniumTNG project that we refer to as ‘MTNG740’ following the naming convention of the IllustrisTNG project. This calculation follows structure formation in a cubic periodic box with side length |$500\, h^{-1}\, \mathrm{Mpc}$| (⁠|$740\, \mathrm{Mpc}$|⁠), using 4320³ dark matter particles of mass |$1.7\times 10^8\, \mathrm{M_\odot }$| and initially 4320³ gas cells with an initial mass of |$3.1\times 10^7\, \mathrm{M_\odot }$|⁠, which is the targeted baryonic mass resolution. This makes MTNG740 15 times larger than TNG300 by volume at a mass resolution that is only 2.8 times worse. The minimum gravitational softening length for gas cells is set to |$\epsilon _\mathrm{gas,min}=370\, \mathrm{pc}$| (it changes with the cell size), while the gravitational softening length of dark matter and stars is set to |$\epsilon _\mathrm{DM,\star }=3.7\, \mathrm{kpc}$|⁠.

MTNG uses the Planck 2016 cosmology (Planck Collaboration XIII 2016a) in identical form as IllustrisTNG in order to facilitate direct comparisons with TNG, i.e. Ω₀ = 0.3089, Ω_nb = 0.0486, Ω_Λ = 0.6911, σ₈ = 0.8159, n_s = 0.9667, and h = 0.6774. We generate the initial conditions at z = 63 with second-order Lagrangian perturbation theory using Gadget4 (Springel et al. 2021). We apply the fixed-and-paired variance suppression technique (Angulo & Pontzen 2016), although we can only afford to run one realization of the pair, which gives however already a good part of the benefit of this approach. For the matching dark matter simulations, we have simulated both pairs, however.

The MTNG740 simulation was executed on 122 880 cores on the SuperMUC-NG machine at the Leibniz Computing Center, for a total wallclock time of 57 d. It consumed 1.7 × 10⁸ core-hours and had a total memory requirement of |$180\, \mathrm{TB}$|⁠. The total data output of the simulation amounts to |$1.1\, \mathrm{PB}$|⁠.

In the following text, we refer to the radius that encloses an average density of 500 times the critical density of the universe as R_500c, and for brevity, we shall usually drop the subscript c in the rest of the paper. We do the same for the analogous quantities M₅₀₀, M₂₀₀, and R₂₀₀.

3 VERIFICATION

Fig. 1 shows a visual impression of the flagship MTNG740 full physics simulation. To verify the galaxy formation results of the MTNG740 full physics box, we compare them to a set of simulations from the IllustrisTNG project. We specifically compare with TNG100, which has the resolution the galaxy formation physics model was calibrated at (Pillepich et al. 2018a), with TNG300, which is the largest box of IllustrisTNG and has a resolution closest to MTNG740, and with TNG300-2, which has the same box size as TNG300 but features an eight times coarser mass resolution. Thus, TNG300 and TNG300-2 bracket MTNG740 in mass resolution. In Table 1, we list the main numerical parameters of these simulations and contrast them with the MTNG740 full physics box. Notably, MTNG740 simulates a nearly 15 times bigger volume than TNG300 at a mass resolution that is worse by a factor of 2.8. This is made possible by following 5.2 times more resolution elements than TNG300.

To assess the galaxy properties computed by the MTNG740 full physics simulation, we reproduce in Fig. 2 the main calibration plots used for the TNG galaxy formation model (Weinberger et al. 2017; Pillepich et al. 2018a) for the four simulations summarized in Table 1. Where applicable, we bin the simulation data logarithmically with four bins per decade. We show data points for individual objects instead of the median of the bin for bins that contain fewer than 10 data points. We show individual objects as points rather than medians once a bin contains fewer than 10 objects.

To guide the eye and to allow a first rough qualitative assessment, we also include updated observational data in every panel, while for the original data used in the calibration of the model we refer to Pillepich et al. (2018a). Note, however, that all data shown here are quantities inferred from observations based on additional assumptions, they are not direct observables. A more powerful comparison would require elaborate forward modelling of mock observations of the simulations, such that synthetic observables could be contrasted directly with observed quantities. This is beyond the scope of this paper but will be addressed in future work.

In the top left panel of Fig. 2, we show the evolution of the cosmic star formation rate density, |$\dot{\rho }_\star \, \mathrm{[M_\odot \, yr^{-1}\, Mpc^{-3}]}$|⁠. We compute it from the initial mass of all star particles formed in the simulations, using 100 logarithmically spaced redshift bins. The star formation rate density of the simulations changes systematically with mass resolution such that simulations with better mass resolution form more stars. MTNG740 lies closest to TNG300, albeit slightly below it at high redshift, and slightly above it at z = 0 where it almost matches TNG100. For comparison, we consider recent results from reconstructions of the star formation rate density of the Universe between z = 0.5 and z = 6 using sub-mm observations (Gruppioni et al. 2020), from infrared observations between z = 0.5 and z = 6 (Rowan-Robinson et al. 2016), and at high redshift from UV data (Oesch et al. 2018). The star formation rate density in MTNG740 is slightly lower than observations for z < 6 but in good agreement with the data points at z = 10. Note that there are known inconsistencies between the star formation rate density reconstructed from star formation rate indicators and the total stellar mass density inferred for the Universe. The latter seems to favour a star formation rate density that is smaller by about a factor of 2 at z = 2 (Wilkins, Trentham & Hopkins 2008; Yu & Wang 2016).

The top centre panel of Fig. 2 shows the galaxy stellar mass function. Here, we measure the stellar mass in a fixed physical radius of |$30\, \mathrm{kpc}$| that is bound to the galaxy according to SUBFIND. Satellite galaxies are included as well. For comparison, we show the galaxy stellar mass function from the Galaxy And Mass Assembly (GAMA) survey (Wright et al. 2017; Driver et al. 2022).

We find good agreement at the high-mass end, |$M_\star \gtrsim 10^{11}\, \mathrm{M_\odot }$|⁠. MTNG740 is lower by about a factor of 2 at fixed stellar mass for galaxies with around |$M_\star \sim 5\times 10^{10}\, \mathrm{M_\odot }$| and about |$30{{\ \rm per\ cent}}$| lower compared to observations for galaxies with stellar masses |$M_\star \lt 10^{10}\, \mathrm{M_\odot }$|⁠.

In the top right panel of Fig. 2, we show the stellar mass–halo mass relation for all central galaxies and their dark matter haloes. We compute the ratio of the total stellar mass within a fixed physical aperture of |$30\, \mathrm{kpc}$| to M₂₀₀ and divide by the cosmic baryon fraction. We compare the simulations to relations inferred from abundance matching that are based on galaxy counts and dark-matter-only simulations of a Λ cold dark matter universe (Moster et al. 2010; Girelli et al. 2020). MTNG740 is in excellent agreement with TNG100 and abundance-matching results for massive haloes with |$M_\mathrm{200}\gt 10^{12.5}\, \mathrm{M_\odot }$|⁠. For less massive haloes with |$M_\mathrm{200}\lt 10^{12.5}\, \mathrm{M_\odot }$|⁠, the stellar mass of galaxies in MTNG740 is smaller at fixed halo mass than for galaxies in TNG100 but comparable to galaxies in TNG300. Note that at these smaller masses there is also large variation in the estimates from abundance matching and other methods that attempt to reconstruct the relation between stellar mass and halo mass from observations (see e.g. Girelli et al.2020).

In the bottom left panel of Fig. 2, we show the relation between the mass of the most massive black hole of a galaxy and its central stellar mass. Following Pillepich et al. (2018a), we measure the stellar mass as the total stellar mass within twice the stellar half-mass radius, R_⋆.1/2. We include central galaxies as well as satellite galaxies. We find good agreement between TNG100, TNG300, and TNG300-2. The black holes in MTNG740 are less massive by up to a factor of 2 compared to the TNG simulations at fixed stellar mass. For comparison, we show the observed relation between black hole mass and stellar mass of the bulge of their host galaxies that was used for calibration of the IllustrisTNG model (Kormendy & Ho 2013), as well as two more recent relations between black hole mass and stellar bulge and core masses (Zhu, Ho & Gao 2021). The new bulge mass relation is essentially identical to the old black hole mass–bulge mass relation. The new core mass relation has the same slope but a slightly higher normalization. We find overall reasonable agreement between MTNG740 and the observed relations. Note, however, that the stellar mass estimate we use here includes a much larger fraction (about |$80{{\ \rm per\ cent}}$|⁠) of the total stellar mass of galaxies than observational estimates of bulge mass. Moreover, at the high-mass end, the stellar mass estimate includes significant contributions from the intercluster light component. We discuss this in a bit more detail in the context of the galaxy size estimates shown in the bottom right panel of Fig. 2. Also note that Borrow et al. (2022) recently identified an unintended behaviour of the black hole repositioning scheme in the code that can cause massive galaxies to lose their central black holes when their gas is completely stripped by a more massive object (e.g. within a galaxy cluster). This is also present in MTNG740, but does not affect the median values shown here.

In the bottom centre panel of Fig. 2, we show the gas fraction within R₅₀₀ and its dependence on M₅₀₀ for haloes in MTNG740, finding good agreement between the different simulations for |$M_\mathrm{200}\gt 5\times 10^{12}\, \mathrm{M_\odot }$|⁠. We compare to observational data for galaxy groups (Lovisari, Reiprich & Schellenberger 2015) and galaxy clusters (Eckert et al. 2016, 2019). We find that MTNG740 is consistent with the latest data points for the most massive clusters, and also agrees on the group scale. For galaxy clusters with masses |$M_\mathrm{200}\gtrsim 5\times 10^{13}\, \mathrm{M_\odot }$|⁠, however, the inferred gas fractions of most observed clusters are significantly below the gas fractions of the MTNG740 clusters. Moreover, the observational data have a much larger scatter at fixed halo mass than the simulations. A more detailed comparison of the gas fractions in TNG300 with observations using mock X-ray images and estimating the gas mass and total mass of galaxy clusters in the more similar way to observations also finds systematically larger values for TNG300 compared to observations, but much better agreement on the scatter at fixed halo mass (Pop et al. 2022). Additional studies will be needed to understand if the differences at |$M_{500}\gtrsim 10^{14}\, \mathrm{M_\odot }$| and in the associated scatter are primarily caused by observational biases, or rather reflect a limitation of the galaxy formation model in the simulations.

Finally, in the bottom right panel of Fig. 2 we show the stellar half-mass radius of galaxies depending on their stellar mass within twice the stellar half-mass radius. The plot includes central as well as satellite galaxies. The stellar half-mass radius is computed for all stars associated with a galaxy by subfind. For the central galaxies of massive haloes this also includes all stars that are part of the halo but not bound to any satellite galaxy, i.e. the intracluster light. This leads to large size estimates for the most massive galaxies. The different simulations are consistent for |$M_\star \gt 10^{11}\, \mathrm{M_\odot }$|⁠, but they show a systematic trend with mass resolution for lower mass galaxies, where better resolution simulations have smaller galaxies at fixed stellar mass. For comparison, we show data from the GAMA survey (Lange et al. 2015) using their radii computed from r-band luminosity profiles, and from 10 000 galaxies of the Manga (Bundy et al. 2015) survey released in SDSS DR17 (Abdurro’uf et al. 2022) using stellar mass and effective radius estimates from the PIPE3D pipeline (Sánchez et al. 2016). The simulations are reasonably consistent with the observations, although galaxies in MTNG740 appear systematically larger than the observed galaxies by about a factor of 2 at fixed stellar mass.

An important further diagnostic of the galaxy population is its colour distribution. In particular, it is fundamental to quantify to which degree our simulated galaxy population reproduces the observed mass-dependent bimodal colour distribution. Since galaxy colour depends strongly on the specific star formation rate this is a good measure of quenching. TNG100 in particular has been shown to agree well with the observed colour distribution (Nelson et al. 2018). In Fig. 3, we compare the distribution of g − r SDSS colours of MTNG740 galaxies to the colour distributions of the TNG simulations. Colours are computed assuming a single stellar population for each star particle with a Chabrier initial mass function (IMF) (Chabrier 2003) and the Bruzual & Charlot (2003) stellar population synthesis model without dust. The luminosity of all star particles bound to a galaxy according to subfind is then summed up to compute the total galaxy luminosity. We find excellent agreement between MTNG740 and TNG100 for low masses |$M_\star \lt 10^{10.5}\, \mathrm{M_\odot }$| and at the highest masses, |$M_\star \gt 10^{11.5}\, \mathrm{M_\odot }$|⁠. The sharp transition that is present in all TNG simulations at a mass scale |$M_\star \sim 5\times 10^{10.5}\, \mathrm{M_\odot }$| is however shifted to higher masses by 0.5 dex in MTNG740, and appears washed out a bit as well. We associate this with differences in the efficiency of the kinetic AGN feedback introduced by the removal of magnetic fields from the TNG physics model (see Pillepich et al. 2018a). The change in colour distribution for galaxies with masses in the range |$10^{10}\, \mathrm{M_\odot }\lt M_\star \lt 10^{10.5}\, \mathrm{M_\odot }$| directly affects the colour-selected galaxy correlation function in this mass range, which agrees less well with observations than TNG300 (Bose et al. 2022).

As a final piece of validation of the MTNG740 model against IllustrisTNG, we look at the baryonic impact on the total matter power spectrum, which is an important prediction of these non-linear simulations and can affect cosmological constraints. In Fig. 4, we show the ratio of the total matter power spectra of the full physics simulations of MTNG740, TNG100, and TNG300 divided by the total matter power spectra of their dark-matter-only counterparts. We find that the baryonic impact of MTNG740 is almost identical to TNG100 and slightly smaller than in TNG300, with a |$1{{\ \rm per\ cent}}$| deviation at |$k=2\, h\, \mathrm{Mpc}^{-1}$|⁠, a minimum just above |$k=10\, h\, \mathrm{Mpc}^{-1}$|⁠, and a relative difference on this scale of |$\sim 20{{\ \rm per\ cent}}$|⁠.

Following van Daalen, McCarthy & Schaye (2020), we find a suppression of the power spectrum at z = 0 from baryonic effects at a scale of |$0.5\, h\, \mathrm{Mpc}^{-1}$| of |$\Delta P/P_\mathrm{DM\, only}=-0.0043$| and an average baryon fraction of all haloes in the range 6 × 10¹³ < M_500c[M_⊙] < 2 × 10¹⁴ at z = 0 of 0.86 in MTNG740. According to the fit in van Daalen et al. (2020), a suppression of |$\Delta P/P_\mathrm{DM\, only}=-0.0034$| for our baryon fraction would be expected, which is roughly consistent but |$25{{\ \rm per\ cent}}$| smaller than our actual suppression value.

We conclude from the results in this section that there is reassuring agreement between MTNG740 and its TNG predecessors, and that most of the results found for the TNG simulations, in particular for TNG100 and TNG300, can also be expected to hold for MTNG740. At some level, one can therefore view it as a very large volume extension of the TNG suite of simulations, just as we intended.

4 GALAXY CLUSTER SCALING RELATIONS

Owing to its large volume, MTNG740 contains a substantial number of massive galaxy clusters. At z = 0 the full physics box encompasses nine galaxy clusters with |$M_{500} \gt 10^{15}\, \mathrm{M_\odot }$|⁠, and more than 2000 galaxy clusters with |$M_{500} \gt 10^{14}\, \mathrm{M_\odot }$|⁠. At z = 0.25, the time of the last full snapshot before z = 0, which is important because observations cover many more clusters at this epoch than at z = 0, the simulation already contains three massive galaxy clusters with |$M_{500} \gt 10^{15}\, \mathrm{M_\odot }$| and more than 1500 galaxy clusters with |$M_{500} \gt 10^{14}\, \mathrm{M_\odot }$|⁠. These numbers allow us to have a detailed look into cluster scaling relations for a representative sample of simulated galaxy clusters, and study their properties in detail.

In Fig. 5, we show six different galaxy cluster scaling relations for all galaxy clusters at z = 0.25 with |$M_\mathrm{500} \gt 10^{14}\, \mathrm{M_\odot }$|⁠. These relations are often used to infer the mass of a galaxy cluster from observables. We compare the galaxy clusters in MTNG740 with a sample of well-studied clusters of the LoCuSS sample (Mulroy et al. 2019) at similar redshifts between z = 0.15 and z = 0.3. This sample provides various scaling relations between cluster observables from multiwavelength observations of 41 galaxy clusters and weak-lensing-based mass estimates. In addition to weak lensing maps, the observations used to derive the scaling relations include X-ray, millimetre, optical, and infrared observations. The weak lensing masses inferred for the sample range from |$2\times 10^{14}\, \mathrm{M_\odot }$| to |$2\times 10^{15}\, \mathrm{M_\odot }$|⁠. Note that even though we use the weak lensing masses as baseline for our comparison, here they also still carry significant systematic uncertainties (Ardila et al. 2021). The mass values M₅₀₀ for the galaxy clusters in MTNG740 are computed based on spherical overdensity measurements around the potential minimum of the galaxy cluster.

The top-left panel of Fig. 5 shows the richness λ. For the MTNG740 galaxy clusters we compute the richness as the number of galaxies within R₅₀₀ of the galaxy cluster with an intrinsic luminosity in the rest-frame r band brighter than |$L_\mathrm{r}\lt L^{*}_\mathrm{r}+1.75=-19.75$|⁠. Here, we assume |$L^{*}_\mathrm{r}=-21.5$| (Montero-Dorta & Prada 2009) to mimic the luminosity cut used to compute the richness for the observed galaxy clusters. At z = 0.25 this roughly corresponds to a stellar mass of |$5\times 10^9\, \mathrm{M_\odot }$| for satellite galaxies in MTNG740, which lies in the regime where the stellar mass of galaxies in MTNG740 is already significantly below the expected stellar mass–halo mass relation, as shown in the top-left panel of Fig. 2. The slope of the richness–mass relation of the MTNG740 clusters agrees reasonably well with observations. The normalization of the richness of the MTNG740 clusters at fixed cluster mass is, however, about a factor of 2 lower than for the observed clusters.

The interpretation of this difference is not obvious, because it could be caused either by a difference in the richness in the MTNG model or the observational measurement, or by a bias in the weak lensing mass estimate of the observed clusters. A reason for a difference in the richness estimate might be the smaller number of galaxies in MTNG740 compared to observations around and below the knee of the galaxy luminosity function, as shown in the top centre panel of Fig. 2. Moreover, not including a model for dust may bias the galaxy luminosities of the MTNG740 cluster galaxies. In contrast to the observational method, our richness estimate depends strongly on the exact value of the luminosity cut. If we alleviate the luminosity cut to L_r < −18 instead, we can make the richness scaling relation match with the observed data points. On the other hand, if the weak lensing masses underestimate the true mass of the observed galaxy clusters by a factor of 2, both relations will also match. However, this would also affect the other scaling relations that agree significantly better (see below) and therefore appear to be a much less likely explanation. Plausibly, a combination of several factors adds up to the discrepancy we see.

In the top centre panel of Fig. 5, we show the Sunyaev–Zeldovich (SZ) scaling relation. We compute the integrated Compton-y parameter in a sphere of radius R₅₀₀ as

where σ_T is the Thomson cross-section, m_e is the electron mass, c is the speed of light, and P_e is the electron pressure, γ = 5/3 is the adiabatic index, X_e = n_e/n_H = 1.158 is the electron-to-hydrogen number density fraction for a hydrogen mass fraction of X = 0.76, and μ = 0.588 is the mean molecular weight for a fully ionized medium with primordial abundances. Here, E_gas is the total thermal energy of the gas within R₅₀₀. For comparison with the observed clusters we show |$D^2_\rm {A} Y_{500}$|⁠, where |$D_\rm {A}$| is the angular diameter distance at z = 0.25 for the MTNG740 cosmology. We find good agreement between MTNG740 and the observed galaxy clusters for both the shape and the normalization of the SZ scaling relation. The MTNG740 clusters show significantly less scatter around the relation than the observed clusters. We postpone a detailed study about the origin of the scatter (Battaglia et al. 2012) and the reason for the deviations from the observational scatter to future work.

The top-right panel of Fig. 5 shows the scaling relation for the core-excised bolometric X-ray luminosity. Detailed X-ray mocks are in principle possible (see e.g. Pop et al. 2022, for TNG300 clusters), but are far beyond the scope of this paper. Here, we compute the bolometric X-ray emissivity of a single cell j following Kannan et al. (2016) and Marinacci et al. (2018) as

where T is the temperature in Kelvin, h is Planck’s constant (not to be confused with the dimensionless value of the Hubble constant used elsewhere in our work), ρ the mass density, and V the volume of a cell (both measured in cgs units). We use |$Z=\sqrt{1.15}$|⁠, g_ff = 1.3, X_i = 1.079, and X_e = 1.16 for a fully ionized gas of primordial composition (Marinacci et al. 2018). We set the integration limits to E_low = 0.7 keV and E_high = 10 keV to match the observations (Mulroy et al. 2019).

Similar to Pop et al. (2022), we exclude cells that have a temperature below |$10^5\, \mathrm{K}$|⁠, star-forming cells, and cells that have a negative net cooling rate (i.e. that are heated). Our estimated X-ray luminosity scaling relation is roughly consistent with observations. The MTNG740 clusters lie on the observed scaling relation, but with significantly smaller scatter. A significant part of the scatter is likely caused by projection effects that become visible when the simulations are better modelled via mock observables (Pop et al. 2022).

In the bottom-left panel of Fig. 5, we show the scaling relation for the total K-band luminosity within R₅₀₀. We compute the total rest-frame K-band luminosity of the MTNG740 clusters as the sum of the individual K-band luminosity of all stars in R₅₀₀. We assume a solar luminosity in the K band of M_{K, ⊙} = 3.39 (Johnson et al. 1966) and do not include any dust model accounting for intrinsic dust in the galaxies and galaxy cluster. The observed data points are computed as the sum of the K-band luminosities of all galaxies of the clusters within a projected radius of R₅₀₀ that are consistent with being a member of the cluster and brighter than an observational detection threshold. The observed K-band luminosities are corrected for Milky Way dust. Both scaling relations agree well with each other in normalization and slope. The scatter of the observed clusters is slightly larger than for the MTNG740 clusters.

The good agreement of the total K-band luminosity indicates that the total stellar mass in R₅₀₀ matches observations well. Since the richness, however, seems to be smaller than observed, this might indicate that galaxies in R₅₀₀ are stripped too quickly in the simulation, so that they contribute the correct stellar mass into R₅₀₀, but at any given time the population of galaxies within R₅₀₀ is fainter than in observations. Note that satellite stripping is known to depend on numerical resolution with lower resolution simulations showing faster stripping of galaxies when falling into more massive haloes (see e.g. Green, van den Bosch & Jiang 2021).

The bottom centre panel of Fig. 5 shows the scaling relation for the K-band luminosity of the central brightest cluster galaxy (BCG). For MTNG740 we use the integrated K-band luminosity of all stars in a fixed physical radius of |$30\, \mathrm{kpc}$| as a proxy for the luminosity of the central BCG. The BCGs in MTNG740 seem to be brighter in the K band than the observed BCGs by about a factor of 2. Note that the scatter and the slope of the scaling relations are roughly consistent with observations. The difference is interesting because it might point to shortcomings of the AGN model and its coupling to the central BCG or excessive mergers of the central galaxy with satellite galaxies. In particular, an inaccurate centring of the black hole has been argued to lead to too massive BCGs (Ragone-Figueroa et al. 2018). Owing to our crude definition of the brightness of the simulated BCGs, we refrain from a detailed investigation of the differences here. We leave this to future work that should be based on a faithful forward modelling of MTNG740, and involve mock images of the simulated clusters including a model for intrinsic dust.

Finally, we show the gas fractions relative to the cosmic baryon fraction in the lower right panel of Fig. 5. For the simulations we compute the gas fraction as the total gas mass in R₅₀₀ divided by total mass in R₅₀₀ (i.e. M₅₀₀). The panel is essentially a zoom-in of the bottom centre panel of Fig. 2 with a different set of data points. The gas fractions of the MTNG740 clusters lie right in the middle of the distribution of the observed data points, but the relation is much tighter for the MTNG740 clusters than for the observed clusters.

For most of the scaling relations we find that the scatter in the observed relations is significantly larger than for the MTNG740 clusters. One potential avenue to explain the differences is observational uncertainties. These contribute to the mass estimate as well as to the observed quantity of the scaling relation, but the error bars usually only include statistical contributions and at best the well-known systematic errors. Moreover, projection effects can play a significant role for the scatter of various cluster observables (Debackere, Hoekstra & Schaye 2022; Pop et al. 2022). Nevertheless, it is also possible that at least some part of this scatter at fixed mass is real and reflects physics not adequately captured in our simulation. In this case, qualitative changes to our physics model will likely be required to significantly increase the object-to-object scatter at fixed mass, for example, for the gas fraction or Compton-y parameter. We emphasize that both better observations and more realistic comparisons between simulations and observations are required to better understand all the sources of the scatter.

5 GALAXY CLUSTER PROFILES

After discussing integrated quantities of galaxy clusters, we now go one step further and focus on their internal radial profiles. We select all 2359 clusters at z = 0 with |$M_\mathrm{500} \gt 10^{14}\, \mathrm{M_\odot }$| and compare them to observed profiles of nearby galaxy clusters from the X-COP project (Eckert et al. 2019; Ettori et al. 2019; Ghirardini et al. 2019). X-COP provides thermodynamic profiles as well as metallicity and stellar mass profiles (Ghizzardi et al. 2021) for 12 nearby galaxy clusters at redshifts 0.04 < z < 0.1, with estimated masses of |$4\times 10^{14}\, \mathrm{M_\odot } \lesssim M_\mathrm{500} \lesssim 10^{15}\, \mathrm{M_\odot }$|⁠. The results are based on deep XMM–Newton X-ray as well as millimetre observations to reconstruct high-precision profiles of the cluster gas. The 12 galaxy clusters of the X-COP survey are selected from the Planck all-sky SZ map and constitute the most significant detections in the corresponding redshift range.

We first look at different thermodynamical profiles of the gas in clusters. Even though these profiles are strongly correlated, they allow us different insights into the thermodynamical state of the gas in the MTNG740 galaxy clusters.

We start with the electron density profile in Fig. 6. We compute the volume-averaged electron density profile in spherical shells and group the MTNG740 clusters into six logarithmic mass bins in the range |$14 \lt \log _{10}(M_\mathrm{500} / \rm {M}_\odot) \lt 15.2$|⁠. We scale the profiles with R₅₀₀ of each cluster, then combine them. We then show the median electron density profile in each mass bin with |$80{{\ \rm per\ cent}}$| percentile bands, and we compare to the observed cluster profiles from the X-COP project (Ghirardini et al. 2019). All cluster profiles compare the 3D profiles of the MTNG740 clusters with the 3D profiles inferred from observed 2D profiles. We find a reassuring agreement between MTNG740 and the X-COP profiles. The typical scatter of the MTNG740 cluster profiles around the median profile is a factor of a few in the centre, and almost zero at radii equal or larger than R₅₀₀. The central electron densities of the median profiles of the MTNG740 clusters tend to be slightly smaller compared to the observed clusters. However, the most extreme clusters of MTNG740 more than cover the range of the observed profiles. The central electron density of MTNG740 clusters increases systematically with cluster mass at fixed scaled radius, but reaches the same value independently of cluster mass in the outer parts, for R ≳ 0.4R₅₀₀, with very small scatter.

We consider pressure profiles in the top panel of Fig. 7, relative to P₅₀₀ defined as

where |$h_{70}=H_0/(70~\rm {km~s}^{-1}\rm {Mpc}^{-1})$| and h(z) = H(z)/H₀, i.e. h(z = 0) = 1 (Nagai, Kravtsov & Vikhlinin 2007). We compute the median pressure profile in the same way as the profile of the electron density, and again compare to the clusters of the X-COP survey (Ghirardini et al. 2019). We find good agreement between the profiles of the MTNG740 clusters and the observed clusters. There is no systematic trend with cluster mass after expressing the pressure in terms of P₅₀₀. The scatter in the profiles of the MTNG740 clusters around R₅₀₀ is slightly smaller than in the observations, but the latter is likely dominated by observational systematics.

In the bottom panel of Fig. 7, we divide the pressure profile by the universal pressure profile (Arnaud et al. 2010) that is often assumed when converting cluster observables to physical properties of the cluster. The electron pressure profiles of the MTNG740 clusters are well described by the universal pressure profile out to R₅₀₀, with deviations smaller than |$20{{\ \rm per\ cent}}$|⁠. However, at larger radii the universal pressure profile systematically underestimates the pressure in the MTNG740 clusters.

We consider the entropy profile in Fig. 8. For the MTNG740 clusters, we define the entropy as |$K=k_\rm {B} T / n_\mathrm{e}^{2/3}$| (Henden et al. 2018) and express it in units of the characteristic entropy |$K_{500} = k_\rm {B} T_{500} / n_\mathrm{e,500}^{2/3}$|⁠, using

and

where |$k_\rm {B}$| is the Boltzmann constant, |$m_\rm {p}$| is the mass of the proton, |$f_\rm {b}=\Omega _\rm {b} / \Omega _\rm {M}$| is the baryon fraction of the Universe, and ρ_c is the critical density of the Universe at z = 0. We again compare to the profiles from the X-COP survey. We find that the median central entropy is significantly lower for more massive clusters in MTNG740. Overall, the entropy profiles of the MTNG740 clusters agree reasonably well with the observed profiles, but the median central entropy of the simulated clusters is systematically higher than the observed central entropy profiles. This can partly be understood as reflecting the differences in the central electron density profile that is lower for the MTNG740 clusters compared to the observed clusters. The kinetic AGN feedback in TNG quickly dissipates and increases the central entropy on the sound crossing time-scale. As a consequence, the cluster population can mostly be described by non-cool-core clusters (with central entropies |$\gt 30\, \rm {keV~cm}^2$|⁠), even though a small number of cool-core clusters exists in MTNG740. Nevertheless, MTNG740 misses out on a large fraction of cool-core clusters (as shown for the TNG300 sample by Barnes et al. 2018). Fitting an entropy profile of the form

in the range 10⁻²R₅₀₀–R₅₀₀ (Cavagnolo et al. 2009; Barnes et al. 2018) we find that 196 or |$8{{\ \rm per\ cent}}$| of all clusters at z = 0 with M₅₀₀ > 10¹⁴M_⊙ have |$K_0\lt 30\rm {keV~cm}^2$| and can be classified as cool-core clusters, consistent with TNG300 (Barnes et al. 2018). None of the nine clusters with M₅₀₀ > 10¹⁵M_⊙ fulfil this criterion. The lack of cool-core clusters is similar to purely thermal AGN feedback models (Altamura et al. 2022), in contrast to light AGN jet feedback that is able to maintain the cool core while self-regulating the cooling ICM (Ehlert et al. 2022; Weinberger et al. 2022).

In Fig. 9, we look at the metallicity profiles of the intracluster gas. For the MTNG740 clusters we first compute the mass-weighted average metallicity in spherical shells for each cluster, then combine the profiles to obtain median profiles in six mass bins. We compare to the observed metallicity profiles by Lovisari & Reiprich (2019) and Ghizzardi et al. (2021). We note that the observations typically measure iron abundance as a proxy for the metallicity of the gas in galaxy clusters. However, we only track the total metal abundance in MTNG740 (forced by memory constraints), and thus use it for the comparison. The median metallicity profiles of MTNG740 clusters for different mass bins are very similar. They show only a small but still systematic trend with the mass of the clusters. Specifically, more massive galaxy clusters have slightly lower gas metallicity at fixed relative radius. This is consistent with basic expectations because the fraction of baryons that is converted to stars and eventually produces metals decreases with halo mass (see also the top right-hand panel of Fig. 2). We find overall good agreement with the observed metallicity profiles. The inner slope of the observed clusters, most notably for the sample of relaxed clusters, seems to be slightly steeper than for the MTNG740 clusters, though a similar trend is not obvious in the X-COP sample. We leave a more detailed analysis that splits the MTNG740 galaxy clusters into relaxed and disturbed clusters in a similar way as done for observed clusters to future work.

Finally, we show cumulative stellar mass profiles in Fig. 10 and again compare to the observed profiles of the clusters of the X-COP project (Ghizzardi et al. 2021). The stellar mass profiles of the MTNG740 clusters computed from their BCG and satellite galaxies look consistent with the stellar mass profiles of the observed clusters that are generated in a similar way. The observed stellar mass profiles are slightly steeper at large radii.

The total stellar mass profiles of the MTNG740 clusters, however, are significantly higher showing that there is a significant amount of stellar mass in a diffuse intracluster light component. This component is missed when only galaxies in the cluster are counted. Interestingly, the total K-band luminosity of the MTNG740 clusters seems to agree well with the scaling relations shown in Section 4 even though the total K-band luminosity of the observed clusters was also calculated from cluster galaxies. We leave this apparent discrepancy for future work.

We find that the median cumulative stellar mass profiles are almost self-similar for different mass bins, with a normalization that increases with cluster mass. There is reasonable agreement of the outer shape of the stellar mass profiles between MTNG740 and observations. However, the inner stellar mass profiles are significantly steeper for the MTNG740 clusters, and their stellar mass is more centrally concentrated than for the observed profiles. The total stellar mass in the cluster is larger in MTNG740 clusters, but the stellar mass in galaxies in MTNG740 clusters is lower than observed. This difference indicates that galaxies are stripped more quickly in MTNG740 than in the real Universe, which can increase the intracluster light contribution in MTNG740 at the expense of the individual galaxies. However, Ahad et al. (2021) argue that, if anything, galaxies are stripped too little in the Hydrangea simulations, which have a roughly similar mass resolution as MTNG740. We compare the properties of the cluster galaxies in MTNG740 with observations in Section 7. Nevertheless, as shown in previous simulations, the intracluster light can contribute a substantial fraction of the stellar mass outside the central BCG (Puchwein et al. 2010). A similar result has later been found in observations with sufficiently deep surface brightness sensitivity (Presotto et al. 2014; de Oliveira, Jiménez-Teja & Dupke 2022). Future work should mimic the stellar light of the simulated MTNG740 clusters and apply the same observational cuts and limitations to facilitate a proper comparison.

6 THE SZ-SIGNAL ON THE LIGHT-CONE

While spherical profiles and integrated quantities on spherical apertures greatly assist physical understanding, these measures cannot be observed directly. Instead, we can only access projections of an object on the sky. More strictly speaking, we only see our past backwards light-cone. To understand the relevance of this difference for SZ observations, we first consider the Compton-y parameter for a full, deep light-cone output of our simulation in Fig. 11. This light-cone extends to redshift z = 5 and covers a square area of 10° × 10° on the sky.

We compute the Compton-y of a pixel on the light-cone in a similar way as Y₅₀₀ in equation (1) from the total thermal energy of the electrons in the volume covered by pixel on the sky as

where E_i is the thermal energy of a cell and D_{A, i} the angular diameter distance at the redshift when the cell crosses the light-cone, Ω is the solid angle of the pixel on the sky, and Σ_i includes all light-cone cells whose centres lie on the pixel.

Massive galaxy clusters clearly stick out with a high central Compton-y up to 10⁻⁴. We find a mean background value on the full light-cone map of |$\bar{y}=1.68\times 10^{-6}$| (mean) and |$\bar{y}=1.12\times 10^{-6}$| (median). This background level is still consistent with the upper bound from Planck after removing galaxy clusters with |$\bar{y}\lt 2.2\times 10^{-6}$| (Khatri & Sunyaev 2015). Note, however, that the latter bound was computed for a slightly different cosmology. Our background level is also consistent with previous estimates from Springel, White & Hernquist (2001), the Magneticum simulation (⁠|$\bar{y}\lt 1.18\times 10^{-6}$|⁠, Dolag et al. 2016), and with estimates from semi-analytical models on top of dark-matter-only simulations (Osato et al. 2018).

In the right panel of Fig. 11, we zoom in on the most massive galaxy cluster at z = 0.25 in the deep light-cone map. This galaxy cluster has a mass of |$M_{500}=7\times 10^{14}\, \mathrm{M_\odot }$| at z = 0.25, and the white circles in Fig. 11 denote its R₅₀₀ and |$5\, R_{500}$| radii. These are typical aperture radii used to estimate the integrated Compton-y parameter of galaxy clusters. Interestingly, we see a great deal of substructure projected on to the cluster, in particular between R₅₀₀ and |$\, 5R_{500}$|⁠, though most of this substructure is in the foreground and background and not in the immediate vicinity of the cluster, as we will show in the following text.

To examine this point further, we consider the differences between measuring the Compton-y on the deep light-cone and from a local aperture that considers the cluster and its large-scale environment extending over 100 Mpc. We compute the projected Compton-y profile of the cluster on the light-cone as shown in the left-hand panel of Fig. 12. This is compared to the projected Compton-y profile of the same cluster when seen in a cylindrical projection of the snapshot time-slice at z = 0.25, using a total depth of |$\sim 100\, \mathrm{Mpc}$| around the cluster centre and the same viewing direction as for the light-cone observer. Both resulting radial profiles are contrasted in the left panel of Fig. 12. Additionally, we also include separately a projected profile of the light-cone, but with contributions only from gas between z = 0.23 and z = 0.26, i.e. to local contributions around the redshift of the cluster with a depth similar to the cylindrical projection.

We find that the profiles computed from the cylindrical aperture or from the local light-cone contributions are essentially identical. In contrast, the radial profile of the full light-cone starts to significantly deviate at a projected radius of ∼R₅₀₀, when the Compton-y approaches the mean background level of the light-cone.

To quantify this difference better we show the integrated Compton-y parameter for the galaxy cluster out to a given radius in the right panel of Fig. 12. We see that Y_{500, Cylinder} is |$10{{\ \rm per\ cent}}$| smaller than Y_{500, Light-cone}, while Y_{500, Sphere} that we used in Fig. 5 is |$20{{\ \rm per\ cent}}$| smaller than Y_{500, Light-cone}.

For a radius of |$5\, R_{500}$| as used by Planck (Planck Collaboration XXVII 2016b) the situation is different. Y_5R500 is essentially identical for the spherical and cylindrical apertures. Thus, the local background that is included in the cylindrical aperture but not the spherical aperture does not significantly contribute to the total signal. However, Y_5R500 measured on the light-cone is 2.5 times larger than measured for either local aperture, as the distant background that can be seen in Fig. 11 contributes very significantly in this case. Hence, a careful background removal as done via matched or scale-adaptive filtering (Tegmark & de Oliveira-Costa 1998; Sanz, Herranz & Martínez-Gónzalez 2001; Herranz et al. 2002a, b; Schäfer et al. 2006; Kay et al. 2012) is crucial. Future studies can use MTNG740 to test the full procedure used by SZ surveys and to identify potential systematics still present in these procedures.

7 GALAXY CLUSTER OBSERVABLES

For all comparisons of galaxy cluster properties with observations discussed so far we have assumed that we precisely know the mass of the observed galaxy clusters. However, mass is unfortunately not a direct observable for galaxy clusters, so mass estimates always require additional assumptions about the galaxy clusters that all introduce significant systematic uncertainties.

We can avoid these assumptions by comparing to observations fully in observational space. Here, we present a first example of this approach with MTNG740 by comparing two directly observable quantities of the MTNG740 clusters with corresponding measurements made for observed galaxy clusters. For definiteness, we chose the richness and the Compton-y parameter. Both quantities are almost completely independent probes of the state and properties of a galaxy cluster. Here, we do not attempt to construct fully synthetic observations that are analysed in the same way as observational data, but we instead compare the theoretical quantities |$Y_{5R_{500}}$|⁠, i.e. the integrated Compton-y parameter in a sphere with a radius |$5\, R_{500}$|⁠, and the richness λ. Both quantities are well defined and easily determinable for the MTNG740 clusters, but they still require some modest assumptions to be made when measuring them from observations. Specifically, we compare the MTNG740 clusters to a sample of galaxy clusters that have both SZ measurements from Planck (Planck Collaboration XXVII 2016b; Planck Collaboration 2017), and a richness measurement from the RedMaPPer catalogue based on the SDSS-8 galaxy catalogue (Rykoff et al. 2014, 2016).

We compare the MTNG740 clusters to the observed sample in Fig. 13 for all MTNG740 clusters at z = 0.25 with minimum mass of |$M_\mathrm{500} \gt 10^{14.4}\, \mathrm{M_\odot }$|⁠, roughly matching the detection limit of galaxy clusters in Planck at this redshift. We compare to the sample of clusters that are detected in Planck and that are part of the RedMaPPer SDSS8 catalogue in the redshift range z = 0.2 and z = 0.3. We compute the richness and integrated Compton-y parameter or the MTNG740 clusters as described in Section 4.

The MTNG740 clusters fall well within the parameter space covered by the observed galaxy clusters. The scatter at fixed richness or fixed Compton-y parameter seems to be larger for the observed clusters than for the MTNG740 clusters. This difference is most easily explained by observational systematics that contribute to the scatter. The MTNG740 clusters lie at the low-richness end of the distribution of observed clusters. However, this offset is much smaller than what we found for the scaling relations (see Fig. 5), which may indicate that the weak lensing mass estimates that we used for the observed clusters could be biased and be responsible for at least parts of the discrepancy.

Importantly, MTNG740 lacks galaxy clusters with the largest values of richness and Compton-y found in the Planck + RedMaPPer sample. A simple explanation for this discrepancy could be that the most massive MTNG740 clusters are significantly less massive than the clusters in the observed sample due to the still limited volume of our simulation. This explanation is consistent with a rough estimate of the observational volume of SDSS-8 that covers roughly |$25{{\ \rm per\ cent}}$| of the sky, which is equivalent to volume about four times larger between z = 0.2 and z = 0.3 than the MTNG740 volume.

To further test this hypothesis, we group both sets of galaxy clusters into four logarithmic mass bins from |$10^{14.4}\, \mathrm{M_\odot }$| to |$10^{15.2}\, \mathrm{M_\odot }$| with a width of |$0.2\, {\rm dex}$| each. We find object counts of the Planck + RedMaPPer sample of 9, 60, 14, and 2, whereas the corresponding counts for the MTNG740 clusters are 191, 72, 19, and 3, for bins of increasing cluster mass. Surprisingly, we find a slightly larger number of objects in the MTNG740 cluster sample at the high-mass end. However, the richness and Compton-y values for MTNG740 clusters are significantly smaller than the values found for Planck + RedMaPPer clusters at a similar estimated mass.

There are two obvious ways to explain the discrepancies. First, if we assume that the SZ-based mass estimate from the Planck SZ catalogue is underestimating the true mass by |$\sim 0.2\, \mathrm{dex}$| we find almost perfect agreement between both samples. In this case the richness, the integrated Compton-y values, as well as the number of clusters in the highest mass bins are consistent within statistical scatter between both sets of clusters after taking into account the differences in probed volumes.

Remarkably, this factor of |$\sim 0.2\, \mathrm{dex}$| to correct the mass estimates of the Planck clusters is very similar to the empirical correction to the SZ mass estimates applied in the ACT SZ catalogue (Hilton et al. 2021). Their correction is based on the observed richness of a subset of the ACT clusters, and a relation between richness and weak lensing mass estimates (McClintock et al. 2019).

One problem with this interpretation is that the MTNG740 clusters do not lie on the richness scaling relation shown in Section 4. This is not completely unexpected as MTNG740 galaxies around and below the knee of the stellar mass function are not massive enough as shown in Fig. 2, even though it is not clear how much this changes the richness estimate. However, it may also indicate more complicated interpretations of the richness estimate of the MTNG740 cluster for different brightness limits of cluster galaxies or possibly an observational bias in the richness estimate of the LoCuSS sample (Mulroy et al. 2019).

As a consistency check of the richness of the MTNG740 clusters, we select all galaxies of all clusters with a mass larger than |$10^{14.6}\, \mathrm{M_\odot }$| that contribute to the richness estimate of Fig. 13 and show their distribution in absolute rest-frame r-band magnitude versus rest-frame SDSS g − r colour space in the top panel of Fig. 14. For the Planck + RedMaPPer sample, we include all SDSS galaxies with a non-zero cluster membership probability and weight them with their individual probability. Moreover, we weight every galaxy with the inverse of the number of clusters in the cluster mass bin of its parent cluster. Here, we use the same four mass bins as in Fig. 13. Note that the SDSS galaxies are de-reddened for Milky Way dust, but not for intrinsic dust. We apply K-corrections (Chilingarian, Melchior & Zolotukhin 2010) to obtain rest-frame r-band magnitudes of the observed galaxies and use the redshift of their parent galaxy cluster to transform them to absolute magnitudes. This yields good agreement in colour and brightness between the MTNG740 cluster galaxies and cluster galaxies of the Planck + RedMaPPer sample, only the scatter in colour is smaller for the MTNG740 cluster galaxies. In the bottom panel of Fig. 14 we show the average number of galaxies per cluster and decade of magnitude for MTNG740 and the Planck + RedMaPPer sample. They are similar at the bright end, but there are progressively fewer galaxies in MTNG740 clusters compared to the Planck clusters the fainter the galaxies become. The lower panel of Fig. 14 also shows that this difference becomes smaller but does not vanish (dashed line) when we assume that the Planck clusters are |$0.2\, \mathrm{dex}$| more massive, which changes which clusters are matched to the MTNG740 clusters. While this is mostly expected from Fig. 13, it could still be explained by statistical variance.

If we alternatively assume that the SZ-based mass estimates of the Planck clusters are correct, we need to conclude that the MTNG740 clusters systematically underestimate the richness as well as the integrated Compton-y. The former is consistent with the richness scaling relation, but the latter would break the SZ scaling relation as shown in Fig. 5. Because Planck measures |$Y_{5R_{500}}$| rather than Y₅₀₀ as shown in the SZ scaling relation, this would require significantly stronger AGN-driven outflows that heat the gas around massive galaxy clusters on scales well beyond R₅₀₀ without changing the properties of the gas within R₅₀₀. Moreover, the Planck + RedMaPPer sample would need to lack a significant number of the most massive clusters expected to be found in the survey volume.

We conclude that comparisons between simulations and observations in observational space are a crucial and promising test of simulations, and they may also help to understand and interpret observational systematics. We thus envision to carry out a full forward modelling of mock light-cone data, and a detailed comparison of the SZ signal of MTNG740 clusters on the light-cone with observations, in future work.

8 SUMMARY AND OUTLOOK

In this introductory paper of the MillenniumTNG project, we have focused on the flagship full physics simulation that evolves a |$500\, h^{-1}\, \mathrm{Mpc}$| (740 Mpc) cosmological box to z = 0 with a baryonic mass resolution of |$3.1\times 10^7\, \mathrm{M_\odot }$| and a physics model close to the IllustrisTNG model. In Section 3, we established that the MTNG740 hydrodynamical simulation is consistent with its IllustrisTNG predecessors as well as with recent observational data (see Fig. 2). In particular, the colour bimodality of galaxies is still present despite the slightly reduced resolution compared to TNG300. The transition from star-forming blue galaxies to quenched red galaxies occurs at a slightly higher stellar mass in MTNG740 than in TNG (see Fig. 3). The baryonic impact on the total matter matter power spectrum in MTNG740 is essentially identical to TNG100 as shown in Fig. 4.

We then presented a first analysis of the integrated properties of galaxy clusters in MTNG740. In Section 4, we discussed different galaxy cluster scaling relations in MTNG740 and compared them to the nearby galaxy cluster sample by Mulroy et al. (2019). We found generally good agreement between MTNG740 and observed clusters. The normalization of the richness of the observed galaxy clusters seems to be about a factor of 2 higher than the richness of the MTNG740 at fixed mass. However, an important caveat is that the mass estimates of the observed clusters are based on weak gravitational lensing mass estimates that might be biased.

We then looked at the internal structure of the MTNG740 clusters in Section 5. We examined internal profiles of hydrodynamical quantities as well as metallicity and stellar mass of the MTNG740 galaxy clusters and compared them to observed profiles of nearby well-observed galaxy clusters from the X-COP survey (Eckert et al. 2019; Ettori et al. 2019; Ghirardini et al. 2019). We again found satisfactory agreement between the profiles of the MTNG740 clusters and observations. We found the main discrepancy in the inner part of the hydrodynamical profiles of clusters that are able to reproduce non-cool-core clusters, but not cool-core clusters, perhaps requiring a more physical AGN feedback model that reproduces both classes of clusters as observed. A possible refinement to the TNG galaxy model could be to include large-scale jets for the most massive systems.

In Section 6, we considered the Compton-y map of a deep light-cone up to redshift z = 5. The average background value of the Compton-y parameter we find is consistent with the latest constraints from Planck. We then compared the Compton-y of a galaxy cluster at z = 0.25 measured in different apertures. We showed that neither spherical nor cylindrical apertures fully describe the signal observed on the light-cone. In particular, for large apertures (⁠|$5\, R_{500}$|⁠) the Compton-y background significantly contributes so that it cannot be ignored and requires matched filtering to derive an unbiased estimate of the cluster-intrinsic Compton-y signal.

Having shown that the MTNG740 galaxy cluster population is overall consistent with observed galaxy clusters but shows some interesting discrepancies, we turned in Section 7 to a comparison that stays fully in observational space and sidesteps the thorny issue of cluster mass estimates. We compared the galaxy clusters in MTNG740 with a set of galaxy clusters found in both, the Planck SZ cluster catalogue, and the RedMaPPer SDSS-8 catalogue. We compared richness and Compton-y, which are both observables that can also be directly extracted from the simulation. Overall, we found quite good agreement in Fig. 13. Moreover, we confirmed in Fig. 14 that the population of cluster galaxies contributing to the richness estimate is consistent between simulation and observations. However, at the same time this has highlighted that at similar richness and Compton-y parameter, the SZ-based M₅₀₀ mass estimates of the Planck cluster catalogue are |$0.2\, {\rm dex}$| higher than the masses measured for the MTNG740 clusters. We discussed possibly explanations for this offset and remark that this discrepancy is completely consistent with the weak lensing and richness-based mass correction applied in more recent SZ cluster surveys (Hilton et al. 2021).

We conclude that MTNG740 offers a great opportunity to study and understand galaxies and galaxy clusters in their large-scale cosmological context. Also, the MTNG740 simulations can be used to interpret data from future large cosmological surveys, and to quantify and potentially correct for observational biases in cosmological measurements. We give further examples of such applications in our set of companion papers, and presently work on additional research in this direction in forthcoming studies.

ACKNOWLEDGEMENTS

We thank the anonymous referee for helpful comments that improved the paper. RP thanks Wilma Trick for help with the colour scheme of the paper. The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time on the GCS Supercomputer SuperMUC-NG at the Leibniz Supercomputing Centre (LRZ) in Garching, Germany, under project pn34mo. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility, with equipment funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1, and ST/S002502/1, Durham University, and STFC operations grant ST/R000832/1. CH-A acknowledges support from the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094-390783311. VS and LH acknowledge support by the Simons Collaboration on ‘Learning the Universe’. LH is supported by NSF grant AST-1815978. CP acknowledges support by the European Research Council under ERC-AdG grant PICOGAL-101019746. SB is supported by the UK Research and Innovation (UKRI) Future Leaders Fellowship (grant number MR/V023381/1).

DATA AVAILABILITY

The MillenniumTNG simulations will be made fully publicly available at https://www.mtng-project.org/ in 2024. The data shown in the figures of this article will be shared upon reasonable request to the corresponding author.

References