bottom right-most panels in Fig. 17, across 50 and 20 kpc, respectively.

Figure 18: Magnetic field strength in the disks (top) and kinematics of TNG50
MW/M31-like galaxies (bottom). In the top panel, the radial profiles of the
mean intensity of the magnetic fields within $\pm 2$ kpc from the mid plane
span from $\lesssim\mu$G levels (particularly in the outer disks or in the
innermost regions of a subset of galaxies) to tens of $\mu$G. Many TNG50
MW/M31 analogs are characterized by B-field strengths that appear compatible
with those inferred for the galaxy (magenta shaded areas) and nearby galaxies
(gray curves, see text for details). In the bottom panel, we quantify the
degree of ordered versus turbulent motion (V${}_{\rm rot}/\sigma$) of the
stellar vs. gaseous component. The gray large circles denote values for each
galaxy obtained considering all its stars, irrespective of age; the smaller
colored circles show the V${}_{\rm rot}/\sigma$ of stars of progressively
younger age (from red to blue): younger stars are in colder orbits, as cold as
those of the H$\alpha$-emitting i.e. star-forming gas.

Much could be quantified about the gaseous disks of MW/M31 analogs, also
thanks to the bounty of the simulated data. We close this overview by
highlighting two aspects of the disk regions of our sample.

Magnetic field strength in the disk region. Starting from vanishingly-small
values at the initial conditions ($10^{-14}$ comoving Gauss at $z>100$,
Pillepich et al., 2018a), magnetic fields in TNG50 galactic disks grow
exponentially at early times owing to a small-scale dynamo effect (Pakmor et
al., 2017) in addition to adiabatic compression and the combined effects of
shear flows and turbulence due to structure formation and stellar and SMBH
feedback processes (Marinacci et al., 2018). As a result, the strength of
magnetic fields in TNG50 MW/M31-like galaxies can be as high as many tens of
$\mu$G, particularly in their innermost regions.

In Fig.18, top, we quantify the $z=0$ magnetic field radial profiles of TNG50
MW/M31-like galaxies in terms of mean strength in cylindrical annuli of $\pm
2$ kpc thickness around the stellar disk plane. Each galaxy is shown with one
curve (blue lines) and all are compared to available inferences of the
magnetic field strength in our Milky Way’s disk (magenta shaded areas, from
Strong et al., 2000; Crutcher et al., 2010), in M101 (dark gray profile
Berkhuijsen et al., 2016), and in five nearby normal NGC galaxies (lighter
gray curves Basu & Roy, 2013). Firstly, although these inferences are
uncertain and the comparison is made at face value, it is reassuring to see
that many TNG50 MW/M31 analogs share similar profiles or B-field values with
observed ones, as is the case for the Auriga simulated galaxies (not shown
here but see Pakmor et al., 2017). Secondly, and differently than in the case
of the Auriga galaxies, not all TNG50 galaxies show declining magnetic fields
strength with distance: many TNG50 systems actually show suppressed magnetic
fields in the inner few kpc from the centers, in tandem with the gas mass
density suppression and inside-out quenching discussed above. In fact, we
showed in previous TNG-based works that stronger magnetic fields organized in
large-scale, ordered, toroidal configurations naturally arise in gas-rich disk
galaxies, i.e. with topologies reflective of the underlying rotationally-
supported gaseous disks (Marinacci et al., 2018) – we quantify the latter
next.

The kinematics of the gaseous and stellar components. Even if the
quantification of e.g. the HI disk heights needs to be postponed to a
dedicated analysis, Fig.18, bottom panel, demonstrates that the gaseous disks
of TNG50 MW/M31 analogs are by far dominated by ordered motions. The maxima of
the rotational velocities exceed the vertical velocity dispersions of the
H$\alpha$-emitting i.e. star-forming gas by factors of many, up to two orders
of magnitude. This is shown with the gray circles (one per galaxy) and with
the values on the x-axis, whereby the V${}_{\rm rot}/\sigma$ are estimated as
in Pillepich et al. (2019).

By contrast, the stellar components of TNG50 MW/M31 analogs are much hotter
than the gas, with V${}_{\rm rot}/\sigma$ ratios about four times lower (gray
circles): these are obtained by measuring the kinematics at $z=0$ of all the
stars in each galaxy, irrespective of age. We had shown that, according to
TNG50, stars are less rotationally-supported than the star-forming gas across
all galaxy stellar masses and epochs (Pillepich et al., 2019), and we
reiterate this here for the case of MW/M31-like galaxies at the present epoch.
In fact, smaller colored circles in Fig.18 show that, in each galaxy and on
average, younger stars are in colder orbits than older stars at the same time
of inspection (blue vs. red): stars born in the last 100 million years are as
rotationally supported, if not more, than the star-forming gas out of which
they form.

It is now understood both observationally and theoretically that, at fixed
stellar mass, galaxies were on average dynamically hotter at earlier epochs
($z\gtrsim 1.5$) than today (see Pillepich et al., 2019, and references
therein). However, Fig.18 alone cannot provide clues as to whether the hotter
components of the stellar disks are such because they were born hotter or
because they underwent heating after birth due to e.g. secular processes or
merger or flyby events. We are fully convinced that both phenomena are, and
have been, at play in TNG50 MW/M31-like galaxies across their evolution: this
is quantified in detail by blueFrankel et al. in prep.

####  3.8.4 Chemical properties of the disks

Whereas stellar radial migration has been often, albeit arguably, evoked as a
possible heating process affecting the stellar kinematics discussed in the
previous Section, its role in determining the metallicity and abundance make
up of the galactic disks is not disputed. In blueBisht et al. in prep we
quantify that, in the average TNG50 MW/M31 analog, disk stars that are
currently located at the Sun-like location have migrated over the last 2
billion years at a rate of 0.5 kpc Gyr-1. In the following, we briefly provide
inputs on the chemical properties of TNG50 MW/M31-like disks, adding to the
zeroth-order characterization of Fig. 13.

The bulk of the TNG50 MW/M31-like galaxies exhibit negative metallicity
gradients, with progressively metal-poorer stars (blueLian et al. 2023) and
gas phases (Hemler et al., 2021) at larger galactocentric distances,
consistently with extra-galactic observations but with large galaxy-to-galaxy
variations. In fact, recently, blueLian et al. (2023) have provided the first
measurement of the global stellar metallicity profile of the Milky Way’s disk,
finding it non monotonic when all stars of different ages are considered: it
has a mildly positive gradient inside a galactocentric radius of 7 kpc and a
steep negative gradient outside. Intriguingly, a few of TNG50 MW-like galaxies
exhibit a similar broken stellar metallicity profile (blueLian et al. 2023).

At any given galactocentric radius, however, the metallicity distribution of
TNG50 MW/M31-like stars can be extremely wide (Chen et al., 2023). Even though
extremely metal-poor stars (e.g. [Fe/H]$<-3$) preferentially occupy the
central regions of MW/M31-like systems, three dozens among the TNG50
MW/M31-like galaxies actually have more than $10^{6.5-7}\,\rm M_{\odot}$ of
extremely metal-poor stars in cold circular orbits (Chen et al., 2023). This
could provide hints to understand the origin of the very metal-poor stars with
near-circular orbits that have been recently observed in the Galaxy (Carollo
et al., 2022, and references therein), also by investigating their various
chemical element abundances. The TNG50 simulated data offers in fact the
possibility of directly following back in time the origin and birth of each
stellar particle in the simulation, in addition to the modeling of stellar
evolution and chemical enrichment: a vast richness of phenomenology awaits to
be uncovered.

###  3.9 The inner regions

In our journey from large to small spatial scales, we have finally reached the
inner regions of TNG50 MW/M31-like galaxies (i.e. $\lesssim$ a few kpc from
the centers).

With their SMBHs, stellar bulges and stellar bars, these centers are the
result too of billions of years of evolution: namely, of mergers with other
galaxies and their contribution of stars, gas and/or SMBHs, of gas inflows
from the disk and from the CGM, of gas outflows triggered by SMBH or stellar
feedback, of stellar radial migration, of dynamical instabilities, etc. They
hence encapsulate enormous information about the overall assembly and history
of their galaxy. For example, according to TNG50, there can be cases among
MW/M31 analogs whereby even the innermost (spherical) 500 pc of the galaxies
are made of up to 20 per cent ex-situ stellar mass, i.e. star accreted from
mergers (Boecker et al., 2022), with another 20 per cent, for the average
MW/M31-like galaxy, composed of stars migrated inward from the outer regions
of the disk (Boecker et al., 2022).

####  3.9.1 On supermassive black holes (SMBHs)

All TNG50 MW/M31-like galaxies host a SMBH at their center.

As it can be seen from the top left panel of Fig.19, the TNG50 simulation
returns a relationship between SMBH mass and galaxy stellar mass that falls in
the ballpark of many, albeit not all, observational constraints for the
general galaxy population (see dotted gray empirical scaling relations). The
SMBHs of TNG50 MW/M31-like galaxies (blue circles) span the range
$10^{7.5-8.7}\,\rm M_{\odot}$. Therefore, whereas the case of Andromeda with
its $\sim 10^{8}\,\rm M_{\odot}$ SMBH (Bender et al., 2005) is well
represented by the TNG50 outcome, no MW/M31 analog hosts a SMBH as small as
that of Milky Way, i.e. $\sim 4\times 10^{6}\,\rm M_{\odot}$ (Bland-Hawthorn &
Gerhard, 2016, and references therein). This is undesirable and reflects a too
tight Magorrian relation in TNG50, in turn possibly due to numerical-
implementation choices for SMBH seeding and feeding that lead to non fully-
realistic outcomes. However, it shall not be considered a dramatically-
limiting failure of the model. Firstly, it is important to point out that,
within the larger-volume but lower-resolution run TNG300, the TNG model does
return many example galaxies of MW-like mass with SMBH masses compatible with
that of SgrA*. Furthermore, on the one hand, even though less unbiased
observational samples seem to return also much smaller SMBHs at fixed galaxy
stellar mass than earlier compilations, to what extent SgrA* of our Galaxy is
an outlier remains a possibility: at face value, it has almost two orders of
magnitude smaller mass than the SMBH of Andromeda, which is only $\sim
0.2-0.3$ dex more massive in stars than the Milky Way. On the other hand, the
impact of SMBHs on galaxies is mostly driven, in the TNG model, by their
feedback and it is hence the energetic of such feedback the essential element
of the SMBH-galaxy interaction. We refer the interested reader to the in-depth
discussion on this issue in Pillepich et al. (2021), Sections 5.3 and 6.5
therein.

Figure 19: The fundamental characteristics of the central regions of TNG50
MW/M31 galaxies. In the top panels, we showcase the SMBH populations of MW/M31
analogs (blue circles), including their mass (left) and their gas accretion
rates, i.e. Eddington ratios (right). We indicate with circles (crosses) the
SMBHs that are in the so-called thermal (kinetic) feedback mode at $z=0$. In
the bottom left, we quantify the bulge-to-total stellar mass ratio of the
sample, contrasted to the stellar mass ratio of pseudo-bulges and connected to
the mass of the central SMBHs (colors). The bottom right panel shows that at
least 2/3 of TNG50 MW/M31-like galaxies exhibit a bar.

In terms of accretion and hence feedback, in fact, the majority of TNG50
MW/M31 analogs host SMBHs that accrete material from the surroundings at very
low rates, typically $10^{-5}-10^{-3}\,\rm M_{\odot}$ yr-1. As it can be
appreciated from Fig.19, right top panel, these SMBHs exhibit Eddington rates
of $\sim 10^{-4}$ on average (see also Fig.11 of Pillepich et al., 2021), and
thus are consistent with the estimates for the Milky Way, with upper limits of
current/past accretion rates of SgrA* at $8\times 10^{-5}\,\rm M_{\odot}$ yr-1
(Quataert et al., 1999).

Within the TNG model, such low Eddington ratios imply that the majority of
SMBHs of MW/M31 analogs exercise feedback in the kinetic mode: 165 of the 198
SMBHs impart high-velocity, intermittent, isotropic, accretion-disk winds (or
small-scale jets: crosses in Fig.19, right top panel). The remaining ones,
mostly at the lower galaxy and SMBH mass end, are in the thermal feedback
mode. The different SMBH feedback in MW/M31-like have been shown to leave
starkly different imprints into the physical properties of the gas in and
around the galaxies (Pillepich et al., 2021; Ramesh et al., 2023), with the
kinetic feedback ultimately being the physical cause for star-formation
quenching (see also previous Sections). Interestingly, however, even though
more than 80 per cent of the TNG50 MW/M31-like galaxies are impacted at $z=0$
by the SMBH kinetic feedback, only about half of them are quenched, as the
feedback effects are not instantaneous. We note that the median SMBH of TNG50
MW/M31-like galaxies transitioned for the first time from thermal to kinetic
feedback mode about $5-6$ billion years ago, but with an enormous galaxy-to-
galaxy variation. Two thirds of the SMBHs that are currently in kinetic mode
have transitioned to thermal mode only once, i.e. have never got back to
sufficiently high accretion rates.

####  3.9.2 On bulges

A detailed analysis of the properties and origin of photometrically-defined
stellar bulges in TNG50 MW/M31-like galaxies has been presented by Gargiulo et
al. (2022), to which the reader is referred. There it was pointed out, among
others, that less than one fifth of the analogs exhibit bulges with high
Sérsic indeces ($n>2$) and that those show, on average, a higher fraction of
ex-situ stars in their kinematically-selected bulges. On the other hand, the
largest majority of TNG50 MW/M31-like galaxies host low Sérsic-index bulges,
and show properties more akin to pseudo-bulges.

We expand here on kinematically-defined bulges of the TNG50 MW/M31 analogs and
plot in Fig. 19, bottom left, the bulge-to-total stellar mass fraction of both
bulges and pseudo-bulges, according to the kinematic decomposition of Zana et
al. (2022). MW/M31 analogs are represented by colored circles, whereas gray
circles denote all TNG50 galaxies within the relevant stellar mass range.

More massive pseudo-bulges seem to cohabitate with more massive classical
bulges. The median bulge-to-total (pseudo bulge-to-total) ratio for TNG50
MW/M31-like galaxies is 0.13 (0.05), well consistent with the constraints for
the Milky Way and Andromeda (magenta and orange shaded regions, from Boardman
et al., 2020b). We find (not shown) that such ratios in TNG50 do not seem to
depend on galaxy stellar mass (at least within the MW/M31-like range). On the
other hand, as the color code demonstrates and as expected, galaxies with more
massive SMBHs tend to exhibit also more massive bulges, a relationship that is
reflected also on that with SFRs – galaxies with low levels of star formation
have higher-than-average bulge-to-total mass ratios (not shown). Suggestively,
the bottom left panel of Fig. 19 also shows that galaxies with more massive
SMBHs, within the MW/M31-mass bin and within the MW/M31-like selection, have
relatively less massive (or even vanishing) pseudo bulges.

On average and for what it pertains bulges, the kinematic decomposition by
Zana et al. (2022) is well consistent with that by Du et al. (2020), according
to which the median bulge-to-total mass fraction for TNG50 MW/M31-like is
about 15, instead of 13, per cent (Sotillo-Ramos et al., 2022). However,
whereas there can be galaxies without pseudo-bulges, the former decomposition
does not return MW/M31-like galaxies without any bulge component. Instead and
because of a different technique to identify kinematically-defined
morphological components, according to the decomposition by Du et al. (2020),
25 of 198 TNG50 MW/M31-like galaxies can be labelled as bulgeless: see Fig. 13
of Sotillo-Ramos et al. (2022). There, we have also shown that a recent rich
merger history does not imply more massive bulges for galaxies that are
nevertheless MW/M31 analogs at $z=0$, contrary to expectations.

Finally, we find no manifest connections between the relative mass of
pseudo/bulges and the presence of stellar bars, discussed next.

Figure 20: Examples of TNG50 MW/M31-like galaxies with a stellar bar. At least
60 per cent of the 198 galaxies in the sample are barred, depending on how the
identification is conducted. Here we show a few examples in face-on and edge-
on projections, in composite stellar light images as in Fig. 2.

####  3.9.3 On stellar bars

At least 60 per cent of the TNG50 MW/M31-like galaxies exhibit a stellar bar
at $z=0$: see Fig. 19, bottom right panel, and a few examples in Fig. 20.

We count 128 TNG50 MW/M31-like galaxies with bars if we adopt the criteria of
Zana et al. (2022) and 150 (among strong and weak ones) if we adopt the
criteria of Rosas-Guevara et al. (2020), of which 118 are in common – we
provide the flags of the bar identification in Section A.2.

Such a high fraction of barred galaxies is consistent with previous results
based on the TNG simulations and with observations, whereby it has been
quantified that galaxies with $M_{*}\sim 10^{11}\,\rm M_{\odot}$ exhibit the
highest relative occurrence of bars at $z=0$ (in comparison to the global
galaxy population, with or without stellar disk). In fact, among disk
galaxies, the fraction of barred galaxies increases with stellar mass both in
observations and in the TNG100 simulation (Rosas-Guevara et al., 2020; Zhao et
al., 2020): among the TNG50 MW/M31 analogs, bars are relatively less frequent
at the MW-mass scale ($\lesssim 10^{10.9}\,\rm M_{\odot}$) than at the
M31-mass scale.

We also find a suggestive connection between bar occurrence and galaxy’s SF
activity, with the fraction of barred TNG50 MW/M31-like galaxies being higher
towards lower levels of SF: Fig. 19, bottom right panel. Moreover and not
surprisingly, the fraction of barred TNG50 MW/M31 analogs is higher among
galaxies with lower gas mass fractions in the disk regions (not shown). This
in line with the TNG100-based results by Rosas-Guevara et al. (2020, 2022) and
Zhou et al. (2020), with the fraction of barred galaxies being higher among
quenched and gas-poor systems. These relationships have been also seen with
observational data, including intriguing correlations between bar occurrence
and past SF activity (Fraser-McKelvie et al., 2020) and hence possibly
interesting connections among bar formation, quenching and SMBH activity (and
mergers: Łokas, 2020).

The bar sizes of TNG50 MW/M31-like galaxies vary between 1 and 10 kpc (also
depending on definition); many bars appear “decoupled” from the spiral arms,
whereas a few exhibit a more typical SBb morphology (i.e. a bar +two-arm
spiral morphology, like NGC1300). Overall, the stellar bars produced by the
TNG50 simulation and their physical properties are reasonably realistic, also
thanks to its enhanced numerical resolution compared to previous simulations:
via a comparison of MaNGA barred galaxies to a selection of barred galaxies
from TNG50, we have shown that TNG50 bars may be shorter but are not slower
than those observed (Frankel et al., 2022), differently than what claimed in
other studies and/or in the case of previous numerical simulations999The
Frankel et al. (2022) results stem from a quantitative and apples-to-apples
comparison between simulated and observed galaxies: firstly, the simulated
galaxies were matched in galaxy stellar mass and size distributions to the
observed ones; furthermore, we have measured bar properties as closely as
possible to what is done with the observed data, including measuring bar sizes
from r-band SDSS-mocked images..

The large sample and the realism of barred MW/M31-like galaxies from TNG50
offer the opportunity to further and conclusively disentangle physical causes
in the bar-quenching-SMBH loop, and to quantify bar formation scenarios and
bar evolution pathways.

##  4 Summary and conclusions

About a decade after the first realistic cosmological hydrodynamical
simulations of single Milky Way-like systems (Fig. 1 and Table 1), TNG50 marks
a new level of maturity and scope for $\Lambda$CDM cosmological simulations of
galaxies.

TNG50 (Nelson et al., 2019b; Pillepich et al., 2019) is a
gravity+magnetohydrodynamics simulation of a large cosmological volume of 50
comoving Mpc per side. It follows the co-evolution of thousands of galaxies
across the mass spectrum, cosmic epoch, and environment (Fig. 6). With about
200 systems that can be considered analogs of our Galaxy and Andromeda (Fig.
2), TNG50 provides a rich and large theoretical dataset with which to learn
about our closest galaxies, and what they tell us about galaxy formation and
evolution in general. Importantly, this is achieved with a physical model that
was not designed nor intended to replicate in detail the properties of our
Galaxy, Andromeda or any other specific observed galaxy, but rather of the
average galaxy population.

In this paper, we have introduced, motivated, and discussed our fiducial
selection for MW/M31 analogs, which we hence propose either for direct usage
or for reference for future analyses. These are identified as those TNG50
galaxies with stellar disky morphology, with a stellar mass in the range of
$M_{*}=10^{10.5-11.2}~{}\rm M_{\odot}$, and within a MW-like Mpc-scale
environment at $z=0$, namely without nearby massive companions and hosted by a
dark-matter halo less massive than $10^{13}\,\rm M_{\odot}$ (Fig. 3 and
Section 2.3).

TNG50 returns 198 of such MW/M31 analogs at $z=0$. These are resolved with
stellar (dark matter) particles of $8.5\times 10^{4}\rm M_{\odot}$ ($4.5\times
10^{5}\rm M_{\odot}$) and with a gravitational softening at low redshift of
about 300 pc. However, thanks to the quasi-lagrangian nature of the underlying
code (AREPO), the spatial resolution in the gas is fully adaptive, with gas
cells as small as $\sim 10-20$ pc (Fig. 5) and computation time steps as short
as 50 kyr in the innermost dense regions of the simulated galaxies. The $\sim
150$ pc average spatial resolution in the star-forming regions (Fig. 5) means
that astrophysical processes such as star formation and stellar feedback are
subgrid only below the scale of giant molecular clouds and that gaseous and
stellar disks as thin as $100-200$ pc are actually realized in TNG50 (Figs.
16, 17 and 20, and Pillepich et al. 2019). At the same time, the hot gaseous
haloes surrounding MW/M31-like galaxies (Fig. 7, Section 3.3.3 and Ramesh et
al. 2023) are filled with numerous small, $\lesssim$ kpc-sized cold clouds
(Nelson et al. (2020) and blueRamesh et al. 2023b). Therefore, thanks to its
unique combination of volume and resolution, TNG50 expands by many factors (2
orders of magnitude) the sample size of cosmologically-simulated MW/M31
analogs with similar ($\times$10 better) numerical resolution (Fig. 1).

Our MW/M31-like sample is deliberately designed to be as complete as possible
– i.e. we did not want to miss any disky galaxies in TNG50 within the target
mass range (Section 2.3.2). However, it may contain a few systems, such as S0
galaxies, that are not in fact analogs of the Galaxy nor Andromeda. In this
paper, we have provided the tools and suggestions on how to make the selection
more stringent, depending on the chosen scientific application.

In fact, any galaxy (real or simulated) looked at in ever-specific and in an
ever-growing number of details will be, to some extent, unique. This is
apparent when contrasting the measured or inferred properties of our Galaxy
and Andromeda (magenta vs. orange shaded areas in many Figures). Even though
our sample selection was designed to accommodate both types of systems
(chiefly via the stellar mass selection, Section 2.3.1) and even though the
two galaxies can be considered similar in the broad context of cluster
galaxies and ultrafaint dwarfs, it is clear that these two galaxies are not
alike in many details. In fact, because our selection is agnostic to the past
assembly history (Fig. 10 and Section 3.4), we can identify many Andromeda-
like galaxies from TNG50. These are, by choice or necessity, rare among zoom-
in simulation efforts to date (see also a discussion on this by Sotillo-Ramos
et al., 2022).

Remarkably, within the TNG50 MW/M31-like sample, it is possible to identify
several simulated galaxies whose many $z=0$ integral and structural properties
are consistent, singly or jointly, with those measured for the Galaxy and/or
Andromeda. The majority are hosted by haloes whose total mass is compatible
with current inferences for the Milky Way (Fig. 8); the sample includes 3
Local Group (LG)-like systems (Section 3.2 and Fig.7); many of them have
satellite populations consistent, in abundance and galaxy properties, with
those of the Galaxy and Andromeda (Section 3.5, Fig. 11, Engler et al. 2021b;
Engler et al. 2022); and a substantial number of TNG50 MW/M31 analogs host one
or more satellites as massive as e.g. the Magellanic Clouds (Section 3.5.1).

The TNG50 MW/M31 analogs are surrounded by gaseous haloes (Section 3.3.3)
whose total CGM mass is consistent with observational estimates for the Milky
Way (Ramesh et al., 2023) and that even exhibit X-ray emitting eROSITA/Fermi-
like bubbles similar in morphology to those detected towards our Galactic
Center (Fig. 9 and Pillepich et al. 2021). Their stellar haloes are also
diverse (Section 3.3.2), with many of them showing 3D power-law slopes of the
stellar mass density profiles consistent with those measured to be very steep
for the Galaxy or more shallow for Andromeda (Sotillo-Ramos et al., 2022).
Moreover, 60 per cent of the TNG50 MW/M31 analogs are star-forming (like the
Milky Way) or just below the star-forming main sequence (like Andromeda; Fig.
12). And their average stellar ages and metallicity are in the ball park of
extra-galactic observations (Fig. 13), albeit varying widely from system to
system.

Even within such a relatively narrow selection, TNG50 reveals a large galaxy-
to-galaxy diversity, not only in terms of global and halo properties, but also
in disk structures (Figs. 14, 17 and Section 3.8) and the innermost regions
(Section 3.9). At least 60 per cent of the TNG50 MW/M31analogs exhibit a bar,
with varying degrees of strength and lengths (Figs. 19, bottom right, and 20).
The total matter surface density at the Sun-like location varies between
20-120 $\rm M_{\odot}$ kpc-2 (16th-84th percentiles; Fig. 15, top) and yet
local dark-matter densities of $0.2-0.5$ GeV cm-3 are a natural outcome of
$10^{10.5-10.9}\rm M_{\odot}$ TNG50 analogs (larger at larger galaxy masses,
and overall consistent with observational estimates; see Fig. 15, bottom).
Magnetic fields in the inner disk regions can be up to a few tens of $\mu$G
strong (Fig.18), and are overall in accordance with estimates for the Galaxy,
but their average strength profiles are not always monotonic.

It is important to emphasize that such a galaxy-to-galaxy diversity is
predicted by TNG50 to be in place also at fixed galaxy stellar mass, i.e. it
is not only the result of a 0.7 dex-wide stellar mass selection. For example,
the stellar disk length (heights) of TNG50 MW/M31-analogs varies by up to a
factor 6 (3) in 0.1 dex bins of stellar mass at the MW-like scale (Fig. 16).
Suggestively, while we do not find a trend between galaxy stellar mass and
bulge-to-total ratio within the sample, galaxies that host more massive SMBHs
are also characterized by more massive bulges, comparable to those of
Andromeda or the Galaxy, respectively (Fig. 19, bottom left).

By combining the outcome of TNG50 with observational estimates taken at face
value, it appears that Andromeda is a rather average disk galaxy for its
stellar mass, whereas the Galaxy lies on the lower envelope of stellar disk
lengths, stellar disk heights, and SMBH masses. TNG50 returns a few galaxies
whose detailed stellar disk structural properties are compatible with the
available measurements for the Galaxy and Andromeda (Section 3.8.2 and Fig.
16), e.g. in stellar disk length (at least as measured so far – see also
Boardman et al. 2020a). However, no TNG50 galaxy exhibits a SMBH mass as small
as that of the Milky Way (Fig. 19, top left).

As this summary implies, in this paper we have deliberately provided a general
and holistic overview of the MW/M31-like galaxies realized by the TNG50
simulation, with the goal of showcasing the diverse number of possible
scientific and astrophysical applications. More in-depth investigations on
specific topics have been highlighted throughout, and many more are ongoing
(see Table 2).

With this paper, we release a series of broadly useful data outputs and
products ($\sim$ 10 TB in total) that build upon the public release of the TNG
simulations (Nelson et al., 2019a). These are described in detail in Section A
and include: the list of IDs and the key properties of our selected TNG50
MW/M31-like galaxies (Section A.2 and Table 4); the list of IDs and main
properties of their galaxy satellites at $z=0$ (Section A.4; and Table 7);
and, importantly, datacubes of 800 comoving kpc per side of all particle-level
simulation data, i.e. of cold dark matter, gas and magnetic fields, stars,
tracer particles, and super massive black holes (Section A.3 and Tables 5 and
6). These are provided for each of the 198 TNG50 MW/M31-like galaxies at $z=0$
and at all available previous points in time along their main progenitor
branch. This public data release aims to maximally facilitate access and
analysis by public users, enabling the community-wide exploitation of the
costly TNG50 simulation for years to come.

## Acknowledgements

AP thanks Neige Frankel, Hans-Walter Rix, Christina Eilers, and Tommaso Zana
for useful inputs; Gergo Popping, Gandhali Joshi, and Nhut Truong for
producing their HI/H2, stellar assembly times, and X-ray emission catalogs,
respectively; and Yetli Rosas-Guevara and Tommaso Zana for making their post-
processing catalogs publicly available. The authors also thank all the
simulator colleagues who have provided the exact number of MW/M31-like
galaxies in their simulations for Table 1. This work was funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID
138713538 – SFB 881 (“The Milky Way System”, subprojects A01 and A06). RR and
DN acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) through
an Emmy Noether Research Group (grant number NE 2441/1-1). RR and DSR are
Fellows of the International Max Planck Research School for Astronomy and
Cosmic Physics at the University of Heidelberg (IMPRS-HD). The TNG50
simulation was run with compute time granted by the Gauss Centre for
Supercomputing (GCS) under Large-Scale Projects GCS-DWAR on the GCS share of
the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart
(HLRS). This analysis has been carried out on the VERA supercomputer of the
Max Planck Institute for Astronomy (MPIA), operated by the Max Planck
Computational Data Facility (MPCDF), where the released data are also stored.

## Data Availability

All the simulation data of the IllustrisTNG runs, including TNG50, are
publicly available and accessible at https://www.tng-project.org/data (Nelson
et al., 2019a). With this paper, we release added-value catalogs and products
for the MW/M31-like galaxies of TNG50, as described in Appendix A: these are
stored and accessible at the permanently available webpage: https://www.tng-
project.org/data/milkyway+andromeda/.

## References

  * Abadi et al. (2003) Abadi M., Navarro J., Steinmetz M., Eke V., 2003, ApJ, 591, 499
  * Agertz et al. (2011) Agertz O., Teyssier R., Moore B., 2011, MNRAS, 410, 1391
  * Agertz et al. (2021) Agertz O., et al., 2021, MNRAS, 503, 5826
  * Anbajagane et al. (2021) Anbajagane D., Evrard A. E., Farahi A., 2021, MNRAS, 509, 3441
  * Anta et al. (2022) Anta P. M. G.-d., et al., 2022, MNRAS
  * Applebaum et al. (2021) Applebaum E., Brooks A. M., Christensen C. R., Munshi F., Quinn T. R., Shen S., Tremmel M., 2021, ApJ, 906, 96
  * Barmby et al. (2006) Barmby P., et al., 2006, ApJ Letter, 650, L45
  * Basu & Roy (2013) Basu A., Roy S., 2013, MNRAS, 433, 1675
  * Bender et al. (2005) Bender R., et al., 2005, ApJ, 631, 280
  * Benisty et al. (2022) Benisty D., Vasiliev E., Evans N. W., Davis A.-C., Hartl O. V., Strigari L. E., 2022, arXiv
  * Berkhuijsen et al. (2016) Berkhuijsen E. M., Urbanik M., Beck R., Han J. L., 2016, Astronomy & Astrophysics, 588, A114
  * Bird et al. (2013) Bird J. C., Kazantzidis S., Weinberg D. H., Guedes J., Callegari S., Mayer L., Madau P., 2013, ApJ, 773, 43
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn J., Gerhard O., 2016, arXiv.org, astro-ph.GA, 529
  * Boardman et al. (2020a) Boardman N., et al., 2020a, MNRAS, 491, 3672
  * Boardman et al. (2020b) Boardman N., et al., 2020b, MNRAS, 498, 4943
  * Boecker et al. (2022) Boecker A., Neumayer N., Pillepich A., Frankel N., Ramesh R., Leaman R., Hernquist L., 2022, MNRAS, 519, 5202
  * Bovy & Rix (2013) Bovy J., Rix H.-W., 2013, ApJ, 779, 115
  * Bovy & Tremaine (2012) Bovy J., Tremaine S., 2012, ApJ, 756, 89
  * Buch et al. (2020) Buch J., Fan J., Leung J. S. C., 2020, Physical Review D, 101, 063026
  * Buck et al. (2020a) Buck T., Obreja A., Maccio A. V., Minchev I., Dutton A. A., Ostriker J. P., 2020a, MNRAS, 491, 3461
  * Buck et al. (2020b) Buck T., Pfrommer C., Pakmor R., Grand R. J. J., Springel V., 2020b, MNRAS, 497, 1712
  * Callingham et al. (2019) Callingham T. M., et al., 2019, MNRAS, 484, 5453
  * Camps & Baes (2015) Camps P., Baes M., 2015, Astronomy and Computing, 9, 20
  * Carollo et al. (2022) Carollo D., Christlieb N., Tissera P. B., Sillero E., 2022, arXiv
  * Ceverino & Klypin (2009) Ceverino D., Klypin A., 2009, ApJ, 695, 292
  * Chemin et al. (2009) Chemin L., Carignan C., Foster T., 2009, ApJ, 705, 1395
  * Chen et al. (2023) Chen L.-H., Pillepich A., Glover S. C. O., Klessen R. S., 2023, MNRAS, 519, 483
  * Chua et al. (2017) Chua K. T. E., Pillepich A., Rodriguez-Gomez V., Vogelsberger M., Bird S., Hernquist L., 2017, MNRAS, 472, 4343
  * Chua et al. (2022) Chua K. T. E., Vogelsberger M., Pillepich A., Hernquist L., 2022, MNRAS, 515, 2681
  * Collins et al. (2011) Collins M., et al., 2011, MNRAS, 413, 1548
  * Costa et al. (2009) Costa E., Méndez R. A., Pedreros M. H., Moyano M., Gallart C., Noël N., Baume G., Carraro G., 2009, The Astronomical Journal, 137, 4339
  * Crain et al. (2015) Crain R., et al., 2015, MNRAS, 450, 1937
  * Crutcher et al. (2010) Crutcher R. M., Wandelt B., Heiles C., Falgarone E., Troland T. H., 2010, ApJ, 725, 466
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C., White S., 1985, ApJ, 292, 371
  * Davé et al. (2016) Davé R., Thompson R., Hopkins P., 2016, MNRAS, 462, 3265
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Deason et al. (2020) Deason A. J., Fattahi A., Frenk C. S., Grand R. J. J., Oman K. A., Garrison-Kimmel S., Simpson C. M., Navarro J. F., 2020, MNRAS, 496, 3929
  * Diemer et al. (2018) Diemer B., et al., 2018, ApJ Supplement Series, 238, 33
  * Donnari et al. (2019) Donnari M., et al., 2019, MNRAS, 485, 4817
  * Donnari et al. (2021a) Donnari M., et al., 2021a, MNRAS, 500, 4004
  * Donnari et al. (2021b) Donnari M., Pillepich A., Nelson D., Marinacci F., Vogelsberger M., Hernquist L., 2021b, MNRAS, 506, 4760
  * Du et al. (2020) Du M., Ho L. C., Debattista V. P., Pillepich A., Nelson D., Zhao D., Hernquist L., 2020, ApJ, 895, 139
  * Dubois et al. (2014) Dubois Y., et al., 2014, MNRAS, 444, 1453
  * Dubois et al. (2021) Dubois Y., et al., 2021, Astronomy & Astrophysics, 651, A109
  * Eilers et al. (2019) Eilers A.-C., Hogg D. W., Rix H.-W., Ness M. K., 2019, ApJ, 871, 120
  * Einasto et al. (1984) Einasto J., Klypin A. A., Saar E., Shandarin S. F., 1984, MNRAS (ISSN 0035-8711), 206, 529
  * Eisert et al. (2022) Eisert L., Pillepich A., Nelson D., Klessen R. S., Huertas-Company M., Rodriguez-Gomez V., 2022, MNRAS, 519, 2199
  * Emami et al. (2021) Emami R., et al., 2021, ApJ, 913, 36
  * Engler et al. (2021a) Engler C., et al., 2021a, MNRAS, 500, 3957
  * Engler et al. (2021b) Engler C., et al., 2021b, MNRAS, 507, 4211
  * Engler et al. (2022) Engler C., Pillepich A., Joshi G. D., Pasquali A., Nelson D., Grebel E. K., 2022, arXiv
  * Fantin et al. (2019) Fantin N. J., et al., 2019, ApJ, 887, 148
  * Fattahi et al. (2016) Fattahi A., et al., 2016, MNRAS, 457, 844
  * Feldmann et al. (2022) Feldmann R., et al., 2022, arXiv
  * Few et al. (2012) Few C. G., Gibson B. K., Courty S., Michel-Dansac L., Brook C. B., Stinson G. S., 2012, Astronomy & Astrophysics, 547, A63
  * Flynn et al. (2006) Flynn C., Holmberg J., Portinari L., Fuchs B., Jahreiß H., 2006, MNRAS, 372, 1149
  * Font et al. (2020) Font A. S., et al., 2020, MNRAS, 498, 1765
  * Frankel et al. (2019) Frankel N., Sanders J., Rix H.-W., Ting Y.-S., Ness M., 2019, ApJ, 884, 99
  * Frankel et al. (2022) Frankel N., et al., 2022, ApJ, 940, 61
  * Fraser-McKelvie et al. (2019) Fraser-McKelvie A., Merrifield M., Aragón-Salamanca A., 2019, MNRAS, 489, 5030
  * Fraser-McKelvie et al. (2020) Fraser-McKelvie A., et al., 2020, MNRAS, 499, 1116
  * Gadotti (2009) Gadotti D. A., 2009, MNRAS, 393, 1531
  * Gallazzi et al. (2005) Gallazzi A., Charlot S., Brinchmann J., White S., Tremonti C., 2005, MNRAS, 362, 41
  * Gallazzi et al. (2021) Gallazzi A. R., Pasquali A., Zibetti S., Barbera F. L., 2021, MNRAS, 502, 4457
  * Garbari et al. (2012) Garbari S., Liu C., Read J. I., Lake G., 2012, MNRAS, 425, 1445
  * Gargiulo et al. (2022) Gargiulo I. D., et al., 2022, MNRAS, 512, 2537
  * Garrison-Kimmel et al. (2017) Garrison-Kimmel S., et al., 2017, MNRAS, 471, 1709
  * Garrison-Kimmel et al. (2019a) Garrison-Kimmel S., et al., 2019a, MNRAS, 487, 1380
  * Garrison-Kimmel et al. (2019b) Garrison-Kimmel S., et al., 2019b, MNRAS, 489, 4574
  * Geehan et al. (2006) Geehan J. J., Fardal M. A., Babul A., Guhathakurta P., 2006, MNRAS, 366, 996
  * Geha et al. (2017) Geha M., et al., 2017, ApJ, 847, 4
  * Genel et al. (2014) Genel S., et al., 2014, MNRAS, 445, 175
  * Gilmore & Reid (1983) Gilmore G., Reid N., 1983, MNRAS, 202, 1025
  * Graczyk et al. (2020) Graczyk D., et al., 2020, ApJ, 904, 13
  * Grand et al. (2017) Grand R., et al., 2017, MNRAS, 467, 179
  * Grand et al. (2021) Grand R. J. J., et al., 2021, MNRAS, 507, 4953
  * Grand et al. (8 12) Grand R. J. J., et al., 2018-12, MNRAS, 481, 1726
  * Guedes et al. (2011) Guedes J., Callegari S., Madau P., Mayer L., 2011, ApJ, 742, 76
  * Guo et al. (2020) Guo R., Liu C., Mao S., Xue X.-X., Long R. J., Zhang L., 2020, MNRAS, 495, 4828
  * Guzmán-Ortega et al. (2023) Guzmán-Ortega A., Rodriguez-Gomez V., Snyder G. F., Chamberlain K., Hernquist L., 2023, MNRAS, 519, 4920
  * Hagen & Helmi (2018) Hagen J. H. J., Helmi A., 2018, Astronomy & Astrophysics, 615, A99
  * Hammer et al. (2007) Hammer F., Puech M., Chemin L., Flores H., Lehnert M., 2007, ApJ, 662, 322
  * Hemler et al. (2021) Hemler Z. S., et al., 2021, MNRAS, 506, 3024
  * Henden et al. (2018) Henden N. A., Puchwein E., Shen S., Sijacki D., 2018, MNRAS, 479, 5385
  * Hirschmann et al. (4 08) Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014-08, MNRAS, 442, 2304
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Huertas-Company et al. (2019) Huertas-Company M., et al., 2019, MNRAS, 489, 1859
  * Joshi et al. (2020) Joshi G. D., Pillepich A., Nelson D., Marinacci F., Springel V., Rodriguez-Gomez V., Vogelsberger M., Hernquist L., 2020, MNRAS, 496, 2673
  * Joshi et al. (2021) Joshi G. D., Pillepich A., Nelson D., Zinger E., Marinacci F., Springel V., Vogelsberger M., Hernquist L., 2021, MNRAS, 508, 1652
  * Jurić et al. (2008) Jurić M., et al., 2008, ApJ, 673, 864
  * Kennicutt (1998) Kennicutt R. C. J., 1998, Annual Review of Astronomy and Astrophysics, 36, 189
  * Khandai et al. (2015) Khandai N., Matteo T. D., Croft R., Wilkins S., Feng Y., Tucker E., DeGraf C., Liu M.-S., 2015, MNRAS, 450, 1349
  * Kraft et al. (2022) Kraft R., et al., 2022, arXiv
  * Leitner & Kravtsov (2011) Leitner S. N., Kravtsov A. V., 2011, ApJ, 734, 48
  * Lelli et al. (2016) Lelli F., McGaugh S. S., Schombert J. M., 2016, The Astronomical Journal, 152, 157
  * Libeskind et al. (2020) Libeskind N. I., et al., 2020, MNRAS, 498, 2968
  * Licquia & Newman (2015) Licquia T. C., Newman J. A., 2015, ApJ, 806, 96
  * Licquia & Newman (2016) Licquia T. C., Newman J. A., 2016, ApJ, 831, 71
  * Lockman & McClure-Griffiths (2016) Lockman F. J., McClure-Griffiths N., 2016, ApJ, 826, 215
  * Lovell et al. (2018) Lovell M. R., et al., 2018, MNRAS, 481, 1950
  * Mao et al. (2021) Mao Y.-Y., Geha M., Wechsler R. H., Weiner B., Tollerud E. J., Nadler E. O., Kallivayalil N., 2021, ApJ, 907, 85
  * Marel et al. (2012) Marel R. P. v. d., Fardal M., Besla G., Beaton R. L., Sohn S. T., Anderson J., Brown T., Guhathakurta P., 2012, ApJ, 753, 8
  * Marinacci et al. (2014) Marinacci F., Pakmor R., Springel V., 2014, MNRAS, 437, 1750
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * Martig et al. (2012) Martig M., Bournaud F., Croton D. J., Dekel A., Teyssier R., 2012, ApJ, 756, 26
  * McConnachie (2012) McConnachie A. W., 2012, The Astronomical Journal, 144, 4
  * McConnachie et al. (2005) McConnachie A. W., Irwin M. J., Ferguson A. M. N., Ibata R. A., Lewis G. F., Tanvir N., 2005, MNRAS, 356, 979
  * Merritt et al. (2020) Merritt A., Pillepich A., Dokkum P. v., Nelson D., Hernquist L., Marinacci F., Vogelsberger M., 2020, MNRAS, 495, 4570
  * Motwani et al. (2022) Motwani B., et al., 2022, ApJ, 926, 139
  * Naiman et al. (2018) Naiman J., et al., 2018, MNRAS, 477, 1206
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, Astrophysical Journal v.462, 462, 563
  * Nelson et al. (2015) Nelson D., et al., 2015, Astronomy and Computing, 13, 12
  * Nelson et al. (2018a) Nelson D., et al., 2018a, MNRAS, 475, 624
  * Nelson et al. (2018b) Nelson D., et al., 2018b, MNRAS, 477, 450
  * Nelson et al. (2019a) Nelson D., et al., 2019a, Computational Astrophysics and Cosmology, 6, 2
  * Nelson et al. (2019b) Nelson D., et al., 2019b, MNRAS, 490, 3234
  * Nelson et al. (2020) Nelson D., et al., 2020, MNRAS, 498, 2391
  * Nelson et al. (2021) Nelson E. J., et al., 2021, MNRAS
  * Nuza et al. (2014) Nuza S. E., Parisi F., Scannapieco C., Richter P., Gottlöber S., Steinmetz M., 2014, MNRAS
  * Pakmor et al. (2017) Pakmor R., et al., 2017, MNRAS, 469, 3185
  * Pietrzyński et al. (2019) Pietrzyński G., et al., 2019, Nature, 567, 200
  * Pillepich et al. (2014) Pillepich A., et al., 2014, MNRAS, 444, 237
  * Pillepich et al. (2015) Pillepich A., Madau P., Mayer L., 2015, ApJ, 799, 184
  * Pillepich et al. (2018a) Pillepich A., et al., 2018a, MNRAS, 473, 4077
  * Pillepich et al. (2018b) Pillepich A., et al., 2018b, MNRAS, 475, 648
  * Pillepich et al. (2019) Pillepich A., et al., 2019, MNRAS, 490, 3196
  * Pillepich et al. (2021) Pillepich A., Nelson D., Truong N., Weinberger R., Martin-Navarro I., Springel V., Faber S. M., Hernquist L., 2021, MNRAS, 508, 4667
  * Planck (2016) Planck 2016, ApJ, 594, A13
  * Popping et al. (2020) Popping G., et al., 2020, ApJ, 891, 135
  * Posti et al. (2019) Posti L., Marasco A., Fraternali F., Famaey B., 2019, Astronomy & Astrophysics, 629, A59
  * Predehl et al. (2020) Predehl P., et al., 2020, Nature, 588, 227
  * Pérez-Montaño et al. (2022) Pérez-Montaño L. E., Rodriguez-Gomez V., Cervantes Sodi B., Zhu Q., Pillepich A., Vogelsberger M., Hernquist L., 2022, MNRAS, 514, 5840
  * Quataert et al. (1999) Quataert E., Narayan R., Reid M. J., 1999, ApJ, 517, L101
  * Ramesh et al. (2023) Ramesh R., Nelson D., Pillepich A., 2023, MNRAS, 518, 5754
  * Riess et al. (2012) Riess A. G., Fliri J., Valls-Gabaud D., 2012, ApJ, 745, 156
  * Rix & Bovy (2013) Rix H.-W., Bovy J., 2013, The Astronomy and Astrophysics Review, 21, 61
  * Robitaille & Whitney (2010) Robitaille T. P., Whitney B. A., 2010, ApJ Letters, 710, L11
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez V., et al., 2015, MNRAS, 449, 49
  * Rodriguez-Gomez et al. (2016) Rodriguez-Gomez V., et al., 2016, MNRAS, 458, 2371
  * Rodriguez-Gomez et al. (2019) Rodriguez-Gomez V., et al., 2019, MNRAS, 483, 4140
  * Rodriguez-Gomez et al. (2022) Rodriguez-Gomez V., et al., 2022, MNRAS, 512, 5978
  * Rosas-Guevara et al. (2020) Rosas-Guevara Y., et al., 2020, MNRAS, 491, 2547
  * Rosas-Guevara et al. (2022) Rosas-Guevara Y., et al., 2022, MNRAS, 512, 5339
  * Salomon et al. (2020) Salomon J.-B., Bienaymé O., Reylé C., Robin A. C., Famaey B., 2020, Astronomy & Astrophysics, 643, A75
  * Samuel et al. (2020) Samuel J., et al., 2020, MNRAS, 491, 1471
  * Sanderson et al. (2020) Sanderson R. E., et al., 2020, ApJ Supplement Series, 246, 6
  * Sawala et al. (2016) Sawala T., et al., 2016, MNRAS, 457, 1931
  * Sawala et al. (2021) Sawala T., McAlpine S., Jasche J., Lavaux G., Jenkins A., Johansson P. H., Frenk C. S., 2021, MNRAS, 509, 1432
  * Scannapieco et al. (2009) Scannapieco C., White S. D. M., Springel V., Tissera P. B., 2009, MNRAS, 396, 696
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Shi et al. (2022) Shi J., et al., 2022, ApJ, 927, 189
  * Sick et al. (2015) Sick J., Courteau S., Cuillandre J.-C., Dalcanton J., Jong R. d., McDonald M., Simard D., Tully R. B., 2015, Galaxy Masses as Constraints of Formation Models, 311, 82
  * Siegel et al. (2002) Siegel M. H., Majewski S. R., Reid I. N., Thompson I. B., 2002, ApJ, 578, 151
  * Sijacki et al. (2015) Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G., Nelson D., Hernquist L., 2015, MNRAS, 452, 575
  * Sivertsson et al. (2018) Sivertsson S., Silverwood H., Read J. I., Bertone G., Steger P., 2018, MNRAS, 478, 1677
  * Smith et al. (2001) Smith R. K., Brickhouse N. S., Liedahl D. A., Raymond J. C., 2001, ApJ, 556, L91
  * Sotillo-Ramos et al. (2022) Sotillo-Ramos D., et al., 2022, MNRAS, 516, 5404
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, Astrophysical Supercomputing using Particle Simulations, 208, 273
  * Springel et al. (2001) Springel V., White S., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Springel et al. (2008) Springel V., et al., 2008, MNRAS, 391, 1685
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Stewart et al. (2008) Stewart K. R., Bullock J. S., Wechsler R. H., Maller A. H., Zentner A. R., 2008, ApJ, 683, 597
  * Stinson et al. (2013) Stinson G., Brook C., Macciò A., Wadsley J., Quinn T., Couchman H., 2013, MNRAS, 428, 129
  * Strong et al. (2000) Strong A. W., Moskalenko I. V., Reimer O., 2000, ApJ, 537, 763
  * Su et al. (2010) Su M., Slatyer T. R., Finkbeiner D. P., 2010, ApJ, 724, 1044
  * Tamm et al. (2012) Tamm A., Tempel E., Tenjes P., Tihhonova O., Tuvikene T., 2012, Astronomy and Astrophysics, 546, A4
  * Tremmel et al. (2017) Tremmel M., Karcher M., Governato F., Volonteri M., Quinn T. R., Pontzen A., Anderson L., Bellovary J., 2017, MNRAS, 470, 1121
  * Treß et al. (2020) Treß R. G., Sormani M. C., Glover S. C. O., Klessen R. S., Battersby C. D., Clark P. C., Hatchfield H. P., Smith R. J., 2020, MNRAS, 499, 4455
  * Truong et al. (2020) Truong N., et al., 2020, MNRAS, 494, 549
  * Truong et al. (2021a) Truong N., Pillepich A., Nelson D., Werner N., Hernquist L., 2021a, MNRAS
  * Truong et al. (2021b) Truong N., Pillepich A., Werner N., 2021b, MNRAS, 501, 2210
  * Varma et al. (2021) Varma S., et al., 2021, MNRAS
  * Vogelsberger et al. (2014a) Vogelsberger M., et al., 2014a, MNRAS, 444, 1518
  * Vogelsberger et al. (2014b) Vogelsberger M., et al., 2014b, Nature, 509, 177
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42
  * Watkins et al. (2010) Watkins L. L., Evans N. W., An J. H., 2010, MNRAS, 406, 264
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, ApJ Supplement Series, 248, 32
  * Wel et al. (2014) Wel A. v. d., et al., 2014, ApJ Letters, 792, L6
  * Wetzel et al. (2016) Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, ApJ Letters, 827, L23
  * Williams (2003) Williams B. F., 2003, The Astronomical Journal, 126, 1312
  * Worthey et al. (2005) Worthey G., España A., MacArthur L. A., Courteau S., 2005, ApJ, 631, 820
  * Zana et al. (2022) Zana T., et al., 2022, MNRAS, 515, 1524
  * Zanisi et al. (2020) Zanisi L., et al., 2020, MNRAS, 501, 4359
  * Zhang et al. (2019) Zhang H., et al., 2019, MNRAS, 484, 5170
  * Zhao et al. (2020) Zhao D., Du M., Ho L. C., Debattista V. P., Shi J., 2020, ApJ, 904, 170
  * Zhou et al. (2020) Zhou Z.-B., Zhu W., Wang Y., Feng L.-L., 2020, ApJ, 895, 92
  * Zhu et al. (2022) Zhu L., et al., 2022, Astronomy & Astrophysics, 660, A20
  * Zinger et al. (2020) Zinger E., et al., 2020, MNRAS, 499, 768
  * Łokas (2020) Łokas E. L., 2020, Astronomy & Astrophysics, 642, L12

##  Appendix A Data products

With this paper, we release a series of value-added and supplementary data
products that build upon the original data release of all the TNG
simulations101010 https://www.tng-project.org/data/ (Nelson et al., 2019a).
With this new material, we aim at providing the community with even richer and
easier-to-use simulation outputs. By design, this release focuses exclusively
on the TNG50 MW/M31-like galaxies selected and studied here.

###  A.1 Preliminary considerations about the existing TNG50 data products and
infrastructure

The TNG50 data is stored and organized exactly as all other TNG simulations
and is already fully publicly accessible: see Nelson et al. (2019b); Pillepich
et al. (2019) and Nelson et al. (2019a), their Section 3, for details. Here we
summarize salient features to orient the novel readers and users.

Snapshots and catalogs. Data associated to the whole TNG50 volume is available
at 100 snapshots, i.e. it has been stored at 100 points in time, between
$z=15$ and $z=0$ (see e.g. Table 3 of Nelson et al., 2019a, or online). The
separation between two subsequent snapshots at low redshift is about 150 Myr.
At any of these points in time, data is provided in a few main flavours:

  * •

particle data, i.e. snapshots;

  * •

halo and galaxy catalogs;

  * •

merger trees;

  * •

supplementary data catalogs.

By particle data we mean properties of all gas cells and of all the DM,
stellar, tracer, and SMBH particles throughout the ($\sim$50 com Mpc)3
simulated volume. These include, for example, their 3D coordinates and
velocities in the reference system of the comoving box, the total
gravitational potential they sample, but also e.g. the internal energy and
magnetic field content in the case of the gas and the mass, metallicity and
metal abundances at birth for the stellar particles, etc. The complete list,
including description and units, is given e.g. in Tables A.4-A.8 of Nelson et
al. (2015), with full updates and descriptions in Nelson et al. (2019a) and
online. For all TNG runs, there are two different types of snapshots: ‘full’
and ‘mini’, which both encompass the entire volume. However, the ‘mini’
snapshots (80 in total) provide only a subset of particle/cell fields
available in the full ones (20 in total).

The halo and galaxy catalogs provide global and integral properties of FoF
haloes and of their central and satellite subhaloes identified by subfind: if
the latter have non-vanishing stellar masses, they are called galaxies. These
catalogs contain quantities like the galaxy stellar mass, the total halo mass,
the instantaneous SFR, and others, on an object-by-object basis. Again, the
full list of provided quantities is given online, at the same 100 points in
time as the snapshots above.

The merger trees allow to connect individually-identified (sub)haloes and
galaxies across cosmic time, i.e. across the snapshots at which data is
stored.

Additional files include technical inputs for I/O, i.e. for accessing subsets
of the snapshot data, which can amount otherwise up to 2.7 TBytes for each
TNG50 snapshot: for example, it is possible to access only the
gravitationally-bound stellar particles of a given galaxy with a given
SubhaloID.

Finally, over the course of the past years, analyses of the simulations data
have allowed the collection and documentation of a plethora of additional
quantities and derived properties, either at the particle or at the
halo/galaxy level. These include, for example, a diversity of measurements of
galaxy sizes and morphological properties based on the particle data or on
photometry mock images, results of kinematic-decomposition algorithms, more
detailed summary statistics of the galaxies’ merger histories, but also the
characterization of the large-scale cosmic web environment. These are all
listed and documented online111111https://www.tng-
project.org/data/docs/specifications/.

Other resolution runs. The flagship TNG50 simulation that is the focus of this
paper is also referred to as TNG50-1, to denote that it represents the
highest-resolution run of the given sampled comoving volume (see Section 2.4).
The same volume (and hence the MW/M31-like galaxies presented here) is also
realized at three lower (i.e. poorer) resolution levels, in the runs called
TNG50-2, TNG50-3, TNG50-4, with mass (spatial) resolution getting 8 (2) times
larger at progressively lower-resolution levels (see Table 1 of Nelson et al.,
2019a).

These additional runs allow us to quantify the effects of numerical resolution
on the scientific and quantitative findings, if needed: see examples of such
studies in the Appendix of Pillepich et al. (2018a, b, 2019); Engler et al.
(2021a) and others. For most scientific applications, we recommend using the
TNG50 MW/M31-like data as provided, i.e. to consider the outcome and the
properties of the TNG50 MW/M31-like galaxies as the best possible realizations
of the TNG model at the TNG50 resolution – this is because model choices and
numerical resolution are to some extent degenerate with each other and not
easy to interpret. In fact, numerical convergence is a complex issue and
working with simulations across multiple resolution levels can be challenging.
However, should numerical-resolution tests be needed, then each TNG50
MW/M31-like galaxy can be contrasted to its lower-resolution counterparts.
Importantly, the SubhaloIDs of the TNG50 MW/M31-like galaxies in TNG50 are not
the same as in the lower-resolution simulations: the mapping is available from
the “Subhalo Matching Between Runs” catalog.

DM-only runs. Moreover, each run of the TNG50 volume has also been performed
including only cold DM and gravity: these are called DM-only runs and allow us
to quantify the effects of baryonic physics and galaxy formation on the phase-
space properies of DM (see e.g. results described in Section 3.3). Also in
this case, the SubhaloIDs of the TNG50 MW/M31-like galaxies between full-
physics and DM-only realizations are not the same, but are provided.

####  A.1.1 TNG50 MW/M31-like galaxies in TNG50 subboxes

For analysis studies that require finer temporal spacing across snapshots than
about $100-200$ Myr, additional and separate cutouts of the full simulated
volume exist for all TNG runs, including TNG50. These are called ‘subboxes’
and are spatial cutouts of fixed comoving location and size.

There are three subboxes in the TNG50 volume, whereby the whole particle data
therein are stored 3600 times between $z=49$ and today, i.e. every few Myr at
recent times. These are visualized also in Fig. 6, as colored shaded cubes in
projection. Subbox0 is meant to include a handful of MW-mass haloes at $z=0$,
with a cubic volume of about 6 comoving Mpc a side; Subbox1 was chosen to
sample a low-density region, with many dwarf galaxies at $z=0$ and a volume of
about 6 comoving Mpc a side; Subbox2 is centered around the most massive TNG50
cluster at $z=0$ and encompasses a cube of about 7.5 comoving Mpc a side. More
details are given in Tables 4 and 5 of Nelson et al. (2019a) and online.

Table 3: Subsample of the 198 TNG50 MW/M31-like galaxies that are located within the boundaries of the high time-resolution subboxes of the TNG50 volume at $z=0$. In the subboxes, the data are stored every few million years (see Tables 4 and 5 of Nelson et al. 2019a), allowing time evolution analyses that require high time output cadence. SubhaloID | Subbox | In Subbox since $z$ | Notes
---|---|---|---
537941 | 0 | 15 | -
543114 | 2 | 0.26 | barred, S0
565089 | 2 | 0.42 | barred
613192 | 2 | 1.36 | barred, S0

Now, as these are comoving spatial cubes out of the whole TNG50 periodic-
boundary volume, individual galaxies may move in and out of these subboxes
along their cosmic evolution. At $z=0$, 4 of the 198 TNG50 MW/M31-like
galaxies are located in subboxes and have been in there since billions of
years: these are listed in Table 3 and are all centrals. An example of time-
evolution study of two of these galaxies exhibiting eROSITA/Fermi-like bubbles
in their CGM can be found in Pillepich et al. (2021).

###  A.2 IDs and properties of TNG50 MW/M31-like galaxies

Table 4: Description of all fields in the “TNG50 MW/M31-like catalog” released with this paper. All fields are float32 unless otherwise specified. The Flag fields are all either 1 or 0. All fields are 198 long: namely, for each field, there is at least one entry per TNG50 MW/M31 analog. All properties refer to $z=0$ and are listed here in conceptual order. Many more properties of the TNG50 MW/M31-like galaxies can be either found in the official Subfind catalogs (by matching the IDs/indexes) or can be computed and measured from the particle-level data also released here (Section A.3). Field |  Units |  Description
---|---|---
SubfindID |  - |  $z=0$ ID of the galaxy, in the space of the Subfind i.e. Subhalo catalog of the whole TNG50 simulation. It is 0 indexed (int32).
FlagCentral |  - |  1 if the galaxy is the central i.e. the first/main/most massive Subhalo object of its FoF halo (int8).
FlagDiskyVisual |  - |  1 if the galaxy appears disky and with spiral arms based on the authors’ visual inspection of 3-band stellar light images, both face-on and edge-on. See Section 2.3.2 (int8).
FlagDiskyMinor2Major |  - |  1 if the stellar-mass minor-to-major axis ratio is smaller than a certain value: $c/a\leq 0.45$. Based on galaxy shapes from Pillepich et al. (2019). See Section 2.3.2 (int8)
FlagDisky |  - |  1 if FlagDiskyVisual=1 OR FlagDiskyMinor2Major=1. Namely 1 if the galaxy has a stellar disky morphology either based on visual inspection or based on the stellar axis ratios. All entries in this catalog should be equal to 1. See Section 2.3.2 (int8)
FlagIsolated |  - |  1 if there is no galaxy with $M_{*}(<30{\rm kpc})\geq 10^{10.5}\,\rm M_{\odot}$ within 500 kpc distance & if $M_{\rm 200c}({\rm host})<10^{13}\,\rm M_{\odot}$. See Section 2.3.3 (int8)
FlagMstars |  - |  1 if $M_{*}(<30{\rm kpc})=10^{10.5-11.2}\,\rm M_{\odot}$. See Section 2.3.1 (int8)
FlagM200c |  - |  1 if $M_{\rm 200c}({\rm host})=6\times 10^{11}-2\times 10^{12}\,\rm M_{\odot}$. See Section 2.3.4 (int8)
FlagNoVirgo |  - |  1 if there is no FoF/Group halo with $M_{\rm 200c}\gtrsim 10^{14}\,\rm M_{\odot}$ within 10 Mpc distance. For TNG50, this includes a constraint on the two most massive systems in the box at $z=0$. See Fig. 6 and Section 3.1 (int8)
FlagS0 |  - |  1 if the galaxy appears as an S0 per visual inspection, namely it has a disky stellar morphology but with no manifest spiral arms. See Section 2.3.2 (int8)
FlagBarredZana22 |  - |  1 if the galaxy has a bar according to Zana et al. (2022). See Section 3.9.3 (int8)
FlagBarredRosas22 |  - |  1 if the galaxy has a bar according to Rosas-Guevara et al. (2020, 2022). See Section 3.9.3 (int8)
FlagMWM31 |  - |  1 if FlagDisky=1 & FlagMstars=1 & FlagIsolated=1. This is the fiducial selection proposed and discussed in this paper. All entries in this catalog should be equal to 1. See Fig. 3. Other more restrictive selections can be imposed using the provided additional flags (int8)
StellarMass_30kpc |  $\rm M_{\odot}$ |  Galaxy stellar mass evaluated by summing up the mass of all stellar particles that are gravitationally bound according to Subfind and within a 3D radius of 30 kpc from the galaxy center. It is the fiducial measure used throughout this paper.
StellarMass_2r1/2 |  $\rm M_{\odot}$ |  Galaxy stellar mass evaluated by summing up the mass of all stellar particles that are gravitationally bound according to Subfind and within a 3D radius of twice the stellar half mass radius. This is equivalent to the 5th entry of SubhaloMassInRadType in the official Subfind catalogs (Nelson et al., 2015, 2019a).
StellarMass_all |  $\rm M_{\odot}$ |  Galaxy stellar mass evaluated by summing up the mass of all stellar particles that are gravitationally bound according to Subfind, with no distance limit. This also includes the mass in e.g. the stellar halo.
HaloMass_M200c |  $\rm M_{\odot}$ |  Total halo mass in terms of the spherical-overdensity measure $M_{\rm 200c}$. For non central galaxies, this is the halo mass of the host.
HaloMass_Mdyn |  $\rm M_{\odot}$ |  Total halo mass obtained by summing up the mass of all particles and cells that are gravitationally bound according to Subfind.
HaloVirialRadius_R200c |  kpc |  Virial radius of the underlying halo in terms of $R_{\rm 200c}$, in analogy with HaloMass_M200c.
SFR_inst |  $\rm M_{\odot}$yr-1 |  Star formation rate of the galaxy based on the instantaneous SFR of the gas, i.e. sum of the individual star formation rates of all gravitationally-bound gas cells in this galaxy. It is equivalent to SubhaloSFR in the official Subfind catalogs (Nelson et al., 2015, 2019a).
SFR_50Myr |  $\rm M_{\odot}$yr-1 |  Time-averaged star formation rate of the galaxy including stars within a 3D aperture of 30 kpc, based on the stellar particles actually produced over the last 50 million years and using their initial mass at birth; as in Pillepich et al. (2019).
DiskScaleLength |  kpc |  Exponential disk length of the galaxy obtained by fitting the stellar mass surface density of stellar particles in disky orbits between 1 and 4 times the stellar half mass radius. See errorbars and fitting procedure in Sotillo-Ramos et al. (2022) and blueSotillo et al. submitted.
DiskScaleHeightThin_8kpc |  pc |  Thin stellar disk height of this galaxy obtained by fitting the vertical stellar mass density distribution of disk stars with a double squared hyperbolic secant functional form in an annulus of $7-9$ kpc from the center. See errorbars, fitting procedure, and additional different measurements in Sotillo-Ramos et al. (2022) and blueSotillo et al. submitted.
DiskScaleHeightThick_8kpc |  pc |  Same as above, for the thick component. Thin and thick disk heights here are meant as geometrical.
SMBH_Mass |  $\rm M_{\odot}$ |  Mass of the central (i.e. most massive) SMBH in this galaxy. The measurement is based on the PartType5/BH_Mass of the official snapshot data (Nelson et al., 2015, 2019a). For MW/M31-like galaxies, this mass is to all effects equivalent to SubhaloBHMass in the official Subfind catalog.
SMBH_AccretionRate |  $\rm M_{\odot}$yr-1 |  Instantaneous gas mass accretion rate into the SMBH of this galaxy, based on SubhaloBHMdot from the official Subfind catalog (Nelson et al., 2015, 2019a).
SMBH_EddingtonRatio |  - |  Instantaneous Eddington ratio of the SMBH of this galaxy.
SMBH_FeedbackMode |  $\rm M_{\odot}$ |  0 if the SMBH of this galaxy is in thermal feedback mode, 1 if it is exercising kinetic feedback mode. See Sections 2.1, 3.6, and 3.9.1 and references therein (int8).
SMBH_CumEnergy_QM |  erg |  Cumulative amount injected into the surrounding gas by the central SMBH in the high accretion-state (quasar) mode, total over its entire lifetime.
SMBH_CumEnergy_RM |  erg |  As SMBH_CumEnergy_QM but for the energy injected in the low accretion-state (wind or kinetic) mode.

With this paper, we release the list of Subhalo i.e. SubfindIDs of the TNG50
MW/M31-like galaxies identified in Section 2.3 and presented throughout this
paper. With this, we invite colleagues to further explore the simulated
galaxies, to apply additional selection criteria as required according to
their targeted scientific application, or to compare to other choices. As
discussed in Section 2.3, the fiducial selection for MW/M31 analogs from TNG50
aimed at returning an as-complete-as-possible sample, even though it may not
be as pure as certain analyses may require.

The table is given in hdf5 format and includes a number of flags and galaxy
properties extracted for this or other analyses, at $z=0$: the list and basic
description are provided in Table 4. Additional integral properties of the
selected galaxies can be obtained by using the particle-level data – see next.

###  A.3 Snapshot data in cutouts for each TNG50 MW/M31 analog

Table 5: Description of the fields in the snapshot cutouts released with this paper (continues in Table 6). There is one hdf5 file for each of the 198 TNG50 MW/M31-like galaxies at each redshift: we provide one cutout per galaxy, each from snapshot 011 to snapshot 099 (i.e. $z=0$) along the main progenitor branch of each galaxy: see text for details. Each cutout contains all particle-level fields present in the original snapshot files across the entire simulated volume (full or mini snapshots depending on the redshift; Section A.1) and contained within a cube centred on the galaxy position (SubhaloPos) and extending $\pm 400$ ckpc along each direction. However, with this release we provide additional quantities for the galaxies under focus that require substantial analysis or post-processing of the whole-simulation data. These additional fields are listed and described below. They are organized by particle type, i.e. PartType = 0, 1, 3, 4, and 5 for gas cells, DM particles, tracers, stars and wind particles, and SMBH particles, respectively. N is the number of each PartType in each file. Field | Dimensions | Units | PartType(s) |  Description
---|---|---|---|---
SubfindID | N | - | all |  The ID of the subhalo that a given particle/cell is a part of at the given redshift, in the space of the Subfind i.e. Subhalo catalog of the whole TNG50 simulation. It is set to -1 if the particle/cell does not belong to any subhalo (int32).
MainSnapshotIndex | N | - | all |  Index into the main snapshot file(s) of this particular particle/cell. This is useful if one wants to load data from one of the post processing quantities that are available online but are not included in these cutouts (int64).
RotatedCoordinates | N,3 | ckpc$/h$ | all |  Spatial position of each particle/cell in the reference system of the main galaxy: a coordinate transformation is first performed to shift to the frame of reference of the (centre of the) main galaxy (SubhaloPos); thus a rotation is performed by diagonalizing the moment of inertia tensor (or mass tensor) such that projecting along the z-axis yields a face-on view of the central galaxy. The mass tensor is computed from stars and star-forming gas within a given aperture: minimum of once (twice) the stellar half mass radius for stars (gas).
RotatedCenterOfMass | N,3 | ckpc$/h$ | PartType0 |  Spatial position of the center of mass of each gas cell, after having applied the translation+rotation transformation as for RotatedCoordinates.
RotatedVelocities | N,3 | km $\sqrt{(}a)/s$ | all |  Spatial velocity of each particle/cell in the reference system where the main galaxy is at rest: a velocity transformation is first performed with respect to the bulk/peculiar velocity of the main galaxy through the simulated volume (SubhaloVel). Then a rotation is applied as for RotatedCoordinates.
RotatedMagneticField | N,3 | see online | PartType0 |  Magnetic field strength of each gas cell after having applied the transformation (rotation only) as for RotatedCoordinates. For the units see description online.
MH | N |  $10^{10}\rm M_{\odot}$/$h$ | PartType0 |  Total neutral hydrogen mass of each gas cell. It equals the sum of the mass of atomic HI plus molecular H2: MH = MHI + MH2 (Popping et al., 2020)
MH2* | N |  $10^{10}\rm M_{\odot}$/$h$ | PartType0 |  Molecular hydrogen mass of each gas cell, according to three different models and hence three different partitioning schemes (denoted in the names with * = ’GK’, ’BR’, ’KMT’), based on Popping et al. (2020); see also Diemer et al. (2018).
Xray_Emission_03_2keV | N | erg cm3/s | PartType0 |  Intrinsic and instantaneous X-ray emission of the gas in the soft broad band [0.3-2] keV. This includes contributions from both the continuum and the lines and is given as X-ray cooling rate. The cooling rate is physically similar to the gas field “GFM_CoolingRate” and it is computed assuming an emission model “APEC” from the XSPEC package (Smith et al., 2001) using the element abundances traced by the simulation (practically it is done by using the VAPEC model in the XSPEC package).Available only at the $z=0$ snapshot. Based on blueTruong et al. 2023 to be submitted.
Xray_Emission_05_2keV | N | erg cm3/s | PartType0 |  Same as above but for the [0.5-2] keV soft X-ray broad band.
Xray_Emission_03_2keV_C | N | erg cm3/s | PartType0 |  Same as Xray_Emission_03_2keV but for continuum emission only.
Xray_Emission_05_2keV_C | N | erg cm3/s | PartType0 |  Same as Xray_Emission_05_2keV but for continuum emission only.
Xray_Emission_Line_X | N | erg cm3/s | PartType0 |  Intrinsic and instantaneous X-ray emission of the gas (as Xray_Emission_03_2keV) in narrow bands at nine specific lines, with X = CV (298.97 eV), CVI (367.47 eV), NVI (419.86 eV), NVII (500.36 eV), OVIIf (560.98 eV), OVIIr (573.95 eV), OVIII (653.49 eV), FeXVII (725.05 eV), and NeX (1021.5 eV). Based on blueTruong et al. 2023 to be submitted.
Table 6: Continues from Table 5, and focuses on the origin and assembly of stellar particles (PartType4). Here the dimension of all fields is equal to N = number of stellar particles in each cutout (we hence omit the column for brevity). Field | Units | PartType(s) |  Description
---|---|---|---
SubfindIDAtFormation | - | PartType4 |  SubfindID of the subhalo i.e. galaxy in which the stellar particle first appeared; $-1$ if it was formed outside of any subhalo (int16).
SnapNumAtFormation | - | PartType4 |  Snapshot number in which the stellar particle first appeared (int16).
InSitu | - | PartType4 |  $1$ if the stellar particle was formed in situ, $0$ if it was formed ex situ, and $-1$ if it does not currently belong to any subhalo. A stellar particle is considered to have been formed in situ if the subhalo in which it was formed lies along the “main branch” of the subhalo in which the stellar particle is currently found. This is decided using the SubLink “baryonic” merger trees (int8). Based on Rodriguez-Gomez et al. 2016.
AfterInfall | - | PartType4 |  1 if the subhalo in which the stellar particle first appeared had already “infalled” into the halo FoF group where it is currently found; 0 otherwise; $-1$ if not applicable i.e., if the particle was formed in situ or if it was formed outside of any subhalo (int16). Based on Rodriguez-Gomez et al. 2016.
AccretionOrigin | - | PartType4 |  This dataset can take the following integer values: 0, 1, and 2 for ex-situ stellar particles that were accreted from completed mergers (i.e., when the subhalo in which the stellar particle formed has already merged with the current subhalo), ongoing mergers (i.e., when the subhalo in which the stellar particle formed has not yet merged with the current subhalo, but will do so at a later snapshot in the simulation), and flybys (i.e., when the subhalo in which the stellar particle formed has not merged with the current subhalo, and will not do so at any future snapshot in the simulation), respectively; and $-1$ if not applicable i.e., if the particle was formed in situ or if it was formed outside of any subhalo (int8). Based on Rodriguez-Gomez et al. 2016.
MergerMassRatio | - | PartType4 |  The stellar mass ratio of the merger in which a given ex-situ stellar particle was accreted (if applicable). The mass ratio is measured at the time when the secondary progenitor reaches its maximum stellar mass. NOTE: this quantity was calculated also in the case of flybys, without a merger actually happening. Based on Rodriguez-Gomez et al. 2016.
DistanceAtFormation | - | PartType4 |  The galactocentric distance when the stellar particle was formed, given in units of the stellar half-mass radius of the parent galaxy at the formation time.
SnapNumAtStripping | - | PartType4 |  The snapshot in which the stellar particle last switched galaxies, and has hence remained in its present galaxy – i.e. time of last stripping. This field is $-1$ for all in-situ stellar particles and for particles that aren’t part of any subhalo (int16).
ProgGalaxyMassAtStripping |  $10^{10}\rm M_{\odot}$/$h$ | PartType4 |  Total stellar mass (as identified by Subfind) of the secondary progenitor galaxy of this stellar particle, from snapshot ’SnapNumAtStripping - 1’.
ProgGalaxyMassInRadAtStripping |  $10^{10}\rm M_{\odot}$/$h$ | PartType4 |  As ProgGalaxyMassAtStripping but for the stellar mass within twice the stellar half mass radius of the secondary progenitor galaxy.
ProgSubhaloMassAtStripping |  $10^{10}\rm M_{\odot}$/$h$ | PartType4 |  As ProgGalaxyMassAtStripping but for the total subhalo mass of the secondary progenitor galaxy, i.e. the sum of all gravitationally-bound mass.
DistanceToHostAtStripping | ckpc$/h$ | PartType4 |  Distance to the centre of the primary galaxy at snapshot SnapNumAtStripping. Here, the centre of the galaxy refers to the position of the most-bound particle of PartType = 4.

For each of the 198 MW/M31-like galaxies at $z=0$, with this paper we provide
a replica of all the TNG50 snapshot data by selecting all particles and cells
within a cube of $\pm$400 comoving kpc from the center of each galaxy. The
latter is denoted as SubhaloPos and is defined as the location of the most
gravitationally bound particle or cell associated to the galaxy, according to
the Subfind algorithm.

We do so for each selected galaxy at both the $z=0$ snapshot as well as across
all previous snapshots up to $z=7$ (snapshot 011), keeping the size of the
cube fixed at $\pm$400 comoving kpc and choosing the center of the box
according to the position of the main branch progenitor of the galaxy – the
main progenitor branch is obtained by using the SublinkGal merger tree
(Rodriguez-Gomez et al., 2015). Redshift $z=7$ is the first snapshot where all
the progenitors of the $z=0$ MW/M31 analogs have a non-zero stellar mass.

Crucially, these cutouts offer a few advantages with respect to the way data
can be accessed with the standard TNG I/O approaches:

  * •

Firstly, we provide all the particle data of, and around, each MW/M31-like
galaxy, namely irrespective of the functioning and results of the Subfind and
FoF halo finders.

  * •

We provide the data of these particular galaxies as a function of time,
without the need to get access to the merger tree files.

  * •

Finally, and most importantly, we augment these cutouts with fields, i.e.
properties, that are either not publicly available or are distributed across
disjoint supplementary catalogs.

In fact, the snapshots of the TNG50 MW/M31 analogs contain all the fields
present in the original snapshot files (full or mini), as mentioned in Section
A.1 and described online.

In addition to the standard fields, we provide essential data of all particles
and cells (coordinates, velocities, magnetic field strength, etc) also in the
coordinate system of the main galaxy, namely we implement a coordinate
translation + rotation such that the center of the cutout is at the center of
the galaxy and the z-axis is perpendicular to the plane of the disk i.e.
aligned with the spin vector of the galaxy. This is measured as the angular
momentum of star-forming gas and stars within twice the stellar half mass
radius of each galaxy.

Additionally, we collect in these cutouts also fields like the post-processed
neutral atomic and molecular hydrogen mass per gas cell, and many quantities
that describe the origin and assembly of stellar particles, such as their in-
situ or ex-situ flags.

All datasets of the snapshot cutouts of the TNG50 MW/M31-like galaxies are
listed and described in Tables 5 and 6.

###  A.4 Properties of the TNG50 MW/M31-like satellites

With this paper and in conjunction with Engler et al. 2022, we also release a
TNG50 MW/M31-like satellite catalog, namely a catalog with the SubhaloIDs and
a selection of properties of all “satellite” galaxies around TNG50 MW/M31-like
galaxies at $z=0$.

This satellite sample comprises all galaxies with $M_{*}\geq 5\times
10^{6}\,\rm M_{\odot}$ and located within a 3D radius of 300 physical kpc from
any of the MW/M31 analog. Objects with SubhaloFlag$\neq$1 are excluded (Nelson
et al., 2019a, and Section A.5). This is the fiducial galaxy selection we have
used in Section 3.5, Fig. 11 and Engler et al. 2021b; Engler et al. 2022 and
comprises a total of 1237 galaxies. It should be noted that this selection
does not depend on the outcome of the FoF and Subfind halo finders and it is
purely based on stellar mass and spatial location. All the particles and cells
of these galaxies are automatically included in the custom data cutouts of
Section A.3.

In Table 7, we describe the provided quantities: these are listed according to
the MW/M31-like galaxy to which they are closest, which is defined as their
host.

Table 7: Description of all fields in the “TNG50 MW/M31-like satellite catalog” released with this paper. All fields are float32 unless otherwise specified. All fields are 1237 long: namely, for each field, there is one entry per satellite of the TNG50 MW/M31 analogs. All properties refer to $z=0$ and are listed here in conceptual order. The same considerations apply as for the TNG50 MW/M31-like main galaxies. Field | Units |  Description
---|---|---
SubfindIDSat | - |  $z=0$ ID of the satellite galaxy, in the space of the Subfind i.e. Subhalo catalog of the whole TNG50 simulation. It is 0 indexed (int32).
SubfindIDHost | - |  $z=0$ ID of the host galaxy of this satellite, in the space of the Subfind i.e. Subhalo catalog of the whole TNG50 simulation. These can only take the values of the SubhaloIDs of the TNG50 MW/M31 analogs of Table 4. The host is the closest main galaxy. It is 0 indexed (int32).
DistanceToHost3D | kpc |  Three-dimensional distance of the satellite galaxy to the closest MW/M31-like hosts.
DistanceToHost2D | kpc |  Two-dimensional distance of the satellite galaxy in a random projection to the closest MW/M31-like hosts (along z axis of the simulation box).
3DVelocityRelToHost | km/s |  Three-dimensional relative velocity of the satellite galaxy with respect to the closest MW/M31-like hosts.
LoSVelocityRelToHost | km/s |  Line-of-sight relative velocity of the satellite galaxy with respect to the closest MW/M31-like hosts, in a random projection (along z axis of the simulation box).
StellarMass_2r1/2 | $\log\rm M_{\odot}$ |  Galaxy stellar mass evaluated by summing up the mass of all stellar particles that are gravitationally bound according to Subfind and within a 3D radius of twice the stellar half mass radius. This is equivalent to the 5th entry of SubhaloMassInRadType in the official Subfind catalogs (Nelson et al., 2015, 2019a).
Magnitude_VBand | mag |  Absolute V-band Buser’s magnitude of the satellite galaxy based on the summed-up luminosities of all the gravitationally-bound stellar particles (Vega magnitudes), from the GFM_StellarPhotometrics field of the official snapshot data (Nelson et al., 2015, 2019a).
Magnitude_rBand | mag |  As for Magnitude_VBand, but in the SDSS r band (AB magnitudes).
Mdyn | $\rm M_{\odot}$ |  Total halo mass of the satellite galaxy obtained by summing up the mass of all particles and cells that are gravitationally bound according to Subfind.
GasMass_2r1/2 | $\log\rm M_{\odot}$ |  Total gas mass of the satellite galaxy evaluated by summing up the mass of all stellar particles that are gravitationally bound according to Subfind and within a 3D radius of twice the stellar half mass radius. This is equivalent to the 1st entry of SubhaloMassInRadType in the official Subfind catalogs (Nelson et al., 2015, 2019a).
HIGasMass_2r1/2 | $\log\rm M_{\odot}$ |  Atomic hydrogen mass of the satellite galaxy evaluated by summing up the HI mass of all gas cells that are gravitationally bound according to Subfind and within a 3D radius of twice the stellar half mass radius. This is based on the ’GK’ HI+H2 partitioning by Popping et al. (2020).
H2GasMass_2r1/2 | $\log\rm M_{\odot}$ |  As HIGasMass_2r1/2 but for the molecular hydrogen.
HalfLightRadius2D_VBand | kpc |  Two-dimensional circularized stellar half-light radius of the satellite galaxy in the V-band, from a face-on projection. Based on Pillepich et al. (2019).
StellarVelDisp_3D | km/s |  Three-dimensional standard deviation of the velocities of all stellar particles of the satellite galaxy within twice the stellar half-mass radius weighted by their respective stellar mass.
SurfaceBrightness2D_rBand | mag arcsec-2 |  Two-dimensional surface brightness of the satellite galaxy in the SDSS r band (see above).
Vmax | km/s |  Maximum of the circular velocity profile of the satellite galaxy, accounting for all matter components.
SFActivity | - |  Flag denoting the star-formation activity of the satellite gaalxy: 1, 0, -1 for star-forming, quenched and green-valley galaxies, respectively. This is based on Pillepich et al. (2019) and on the distance from the ridge of the star-forming main sequence (int8).
t10 | Gyr ago |  Stellar assembly time $\tau_{10}$, i.e. time when the galaxy assembled 10 per cent of its $z=0$ stellar mass. Based on Joshi et al. (2021).
t50 | Gyr ago |  As t10 but for the time when the galaxy assembled 50 per cent of its $z=0$ stellar mass.
t90 | Gyr ago |  As t10 but for the time when the galaxy assembled 90 per cent of its $z=0$ stellar mass.

###  A.5 Scientific remarks and other practical inputs

Any user of the data that is released with this paper must be aware of the
possible limitations and areas of caution that have been identified in past
analyses. In fact, the TNG50 simulation, and the TNG model in general, have
been shown to return realistic outcomes in many astrophysical regimes and to
various levels of qualitative or quantitative agreement (see Sections 1 and
2.1 and references therein). The results put forward and referenced to in this
paper add to this encouraging evidence for the particular case of MW/M31-like
galaxies: in particular, they demonstrate that TNG50 is capable of returning
systems whose many integral and structural properties are consistent, one or a
few at the time, with those measured for two particular galaxies, the Milky
Way and Andromeda. However, numerous areas still require additional study: the
comparison with observational results is inherently complex, and the
underlying galaxy formation model is necessarily simplified and hence involves
limitations and uncertainties.

In the TNG data release paper (Nelson et al., 2019a), we have identified
regimes of caution for the usage of the data and the interpretation of the
ensuing scientific results: see their Section 5. Below we stress upon already-
known tensions and limitations and include additional considerations relevant
to the scientific exploration of MW/M31-like analogs with TNG50:

  * •

Cold and dense gas. In the TNG simulations, and hence in TNG50, gas colder
than $10^{4}$ K is not modeled and hence not realized. Gas denser than the
effective threshold for star formation (0.13 atoms cm-3) is labeled as star
forming and its temperature is ‘effective’ i.e. intermediate between cold and
hot phases (see Section 2.1 and references therein). Studies with TNG50 of the
cold and dense ISM are therefore limited. Forms of molecular hydrogen and the
atomic-to-molecular transitions can be however modeled in post processing (see
e.g. Diemer et al., 2018; Popping et al., 2020) and we have included those
estimates in this release (see Table 5).

  * •

Stellar particles. Stellar particles in TNG50 represent mono-age stellar
populations of about $8\times 10^{4}\,\rm M_{\odot}$: although the model
accounts for their aging, mass return and metal production and return, stellar
particles do not represent individual stars. Mocks of stellar catalogs
targeting specific surveys, such as with Gaia, can be however constructed in
post processing (see e.g. Grand et al., 8 12; Sanderson et al., 2020, for
other simulations).

  * •

Subhaloes of non-cosmological origin. As extensively discussed by Nelson et
al. (2019a, Section 5.2 therein), not all objects identified by Subfind and
stored in the corresponding catalogs should be considered ’galaxies’. This is
particular relevant when considering satellite galaxies, as some of them may
not have formed in their own DM halo and have not a “cosmological” origin.
Such objects can be easily discarded with the SubhaloFlag provided in the
official TNG data release and in the cutouts of Section A.3. In stellar and
gas surface mass density maps and stellar light maps, these objects appear
like massive star clusters (e.g. in Fig. 14). However, it is not clear whether
their formation is physical or due to numerical artifacts and so users should
treat them with caution. A quantitative analysis of their occurrence can be
found in Boecker et al. (2022).

  * •

UV background at high z. In all the TNG simulations, and hence also in TNG50,
the UV background is switched on at $z=6$ instead of following the gradual
build-up since higher redshifts. This may impact the progenitors of TNG50
MW/M31-like galaxies and their satellites during the epoch of Reionization.

  * •

Alpha abundances. As shown and discussed by Naiman et al. (2018), the
simulated stars in MW-like galaxies of the TNG100 simulation are too alpha-
enhanced (specifically for [Mg/Fe]) in comparison to observations of our
Galaxy. This is the case also in TNG50 and may be due to either too-low SNIa
rates or too high $\alpha$-element yields from SNII adopted in the fiducial
TNG model. As SNIa only contribute to the metal abundances, as yield stellar
tables are notoriously uncertain, and as metals in the TNG model are passive
scalars, i.e. they do not dynamically affect the physical evolution of the
systems, the mismatch does not impact the functioning and outcome of the TNG50
simulation in general. In fact, analyses relying or focusing on abundance
ratios can still be worthy and useful, especially when studied in relative
terms and keeping in mind that they depend on the underlying choices of the
yield tables (see Pillepich et al., 2018a, and Table 2 and Fig. 1 therein). A
global rescaling in post processing of the alpha abundances is a possible way
forward, especially given that the needed modifications to bring TNG50 stars
in agreement with the measurements of the Milky Way are all acceptable within
current observational constraints for SNIa rates and stellar yield tables
(Naiman et al., 2018).

  * •

SMBH masses. As shown and discussed in Section 3.9.1, there are no galaxies
within the TNG50 MW/M31 analogs whose SMBH have a mass similar to that of
SgrA*: scientific applications that rely on the specific mass value of the
central SMBH may be limited.

To streamline and guide the usage of the data released here, we also note the
following points121212Many additional practical inputs can be obtained also in
the Discussion Forum and FAQs at https://www.tng-project.org/data/.:

  * •

Units. The quantities in the data cutouts (Section A.3, Tables5 and 6) are
given in code/simulation units for consistency with the particle-level fields
present and documented in the original snapshot files: useful notes on
conversions can be found online.

  * •

Metal abundances. Individual abundances of nine species (H, He, C, N, O, Ne,
Mg, Si, Fe) are provided for both gas cells and stellar particles. These are
all metal mass fractions, i.e. the dimensionless ratios of mass in the given
species to the total gas-cell or stellar-particle mass (see online
documentation). An additional conversion is required to obtain number
fractions in solar units, e.g. [Fe/H] or [Mg/Fe], where for example [Fe/H] =
$\mathrm{log}_{10}(N_{\mathrm{Fe}}/N_{\mathrm{H}})-\mathrm{log}_{10}(N_{\mathrm{Fe,\odot}}/N_{\mathrm{H,\odot}})$,
where $N_{\mathrm{Fe}}$ and $N_{\mathrm{H}}$ are the fractional abundances of
iron and hydrogen, respectively, and $N_{\mathrm{Fe,\odot}}$ and
$N_{\mathrm{H,\odot}}$ correspond to the solar abundances.

  * •

Metal tagging. The production and mixing of the sum of all metals produced by
AGB stars, SNIa and SNII are separately followed in the TNG simulations and
recorded for both gas cells and stellar particles (Pillepich et al., 2018a).
Additionally, we keep track of the iron mass separately “produced” by SNIa and
SNII (see e.g. maps in Fig. 14). These are convenient features of the released
data but should be used and interpreted with caution, as these tags do not
track where a given heavy element was created, but rather identify the last
star it was ejected from (see notes in Nelson et al., 2019a).

  * •

Europium. The TNG model also includes the injection of r-process material from
neutron stars–neutron star mergers, which are modeled similarly as SNIa
(Pillepich et al., 2018a; Naiman et al., 2018). The total metal mass produced
by NSNS is also tracked and stored for both gas cells and stellar particles.
This allows to extract estimates of e.g. Europium: see Naiman et al. (2018)
and online documentation.

  * •

Stellar quantities at birth and stellar ages. To circumvent the limitations of
the time-wise sparsity of the data outputs, a few quantities are stored and
provided for each stellar particle, such as the spatial position and velocity
at the time of birth, particularly relevant to study e.g. radial migration.
The latter is also provided so that the ages i.e. the formation times of the
stellar particles are available irrespective of the time cadence of the output
data. The temporal resolution of the stellar ages in TNG50 is as good as $\sim
0.03$ Myr, and permits to extract star formation histories with exquisite time
resolution (as in Figs. 12 and 10).

  * •

X-ray emission quantities. The novel and additional X-ray emission quantities
provided in the cutouts (Table 5) are given as cooling rates. Here is an
example of how to compute, e.g., the X-ray luminosity of each gas cell in the
[0.5-2] keV band: $L_{x}[\rm{erg~{}s}^{-1}]={\rm
cooling~{}rate}[\rm{erg~{}cm}^{3}$ s${}^{-1}]\times n_{e}[cm^{-3}]\times
n_{h}[$cm${}^{-3}]\times V[$cm${}^{3}]$, where $n_{e}$ and $n_{h}$ are
hydrogen and electron number densities, and V is the volume of the gas cell.

  * •

Shock finder in subboxes. Another attractive code feature of the TNG
simulations, including TNG50, is the on-the-fly shock finder, which returns
the energy dissipated via shocks and their Mach number. However, although the
two shock finder fields exist also in the subboxes (Section A.1.1), they are
not updated for each subbox output, but only occasionally: hence these two
fields in the subboxes should not be used at present.

  * •

Additional useful catalogs. A number of post-processing catalogs exist for the
TNG50 simulation in general: these include additional quantities related to
the TNG50 MW/M31-like sample of this paper. For example: various
kinematic/morphological decompositions, SDSS-like spectra, MaNGA-like
intregral-field-unit mocks, photometry for tens of broadbands from FUV to
submillimetre, etc. These are available online.