The merger and assembly histories of Milky Way- and M31-like galaxies with TNG50: disk survival through mergers

We analyze the merger and assembly histories of Milky Way (MW) and Andromeda (M31)-like galaxies to quantify how, and how often, disk galaxies of this mass can survive recent major mergers (stellar mass ratio $\ge$ 1:4). For this, we use the cosmological magneto-hydrodynamical simulation TNG50 and identify 198 analog galaxies, selected based on their $z=0$ stellar mass ($10^{10.5-11.2} {\rm M_{\odot}}$), disky stellar morphology and local environment. Firstly, major mergers are common: 85 per cent (168) of MW/M31-like galaxies in TNG50 have undergone at least one major merger across their lifetime. In fact, 31 galaxies (16 per cent) have undergone a recent major merger, i.e. in the last 5 Gyr. The gas available during the merger suffices to either induce starbursts at pericentric passages or to sustain prolonged star formation after coalescence: in roughly half of the cases, the pre-existing stellar disk is destroyed because of the merger but reforms thanks to star formation. Moreover, higher merger mass ratios are more likely to destroy the stellar disks. In comparison to those with more ancient massive mergers, MW/M31-like galaxies with recent major mergers have, on average, somewhat thicker stellar disks, more massive and somewhat shallower stellar haloes, larger stellar ex-situ mass fractions, but similarly massive kinematically-defined bulges. All this is qualitatively consistent with the different observed properties of the Galaxy and Andromeda and with the constraints on their most recent major mergers, 8-11 and ~2 Gyr ago, respectively. According to contemporary cosmological simulations, a recent quiet merger history is not a pre-requisite for obtaining a relatively-thin stellar disk at $z=0$.


INTRODUCTION
The Lambda cold dark matter (ΛCDM) cosmological model successfully explains the evolution of the Universe and the formation of the large-scale cosmic web. Structures grow hierarchically, with more massive haloes growing by smooth accretion and by merging with smaller ones (Genel et al. 2010), in a process dominated by (cold) dark matter, down to galactic scales.
In this context, mergers and interactions between galaxies are expected to be a frequent phenomenon, more so the higher the redshift (Fakhouri & Ma 2008). In fact, galaxy mergers are known to play a central role in the mass growth of galaxies. Within the hierarchical growth of structure scenario, more massive haloes and galaxies undergo, on average, a larger number of mergers than lower-mass ones (Fakhouri et al. 2010;Rodriguez-Gomez et al. 2015). Moreover, more massive galaxies have been shown, via cosmological hydrodynamical galaxy simulations, to be made of larger fractions of stellar material that is accreted from mergers and orbiting satellites (e.g. Rodriguez-Gomez et al. 2016;Pillepich et al. 2018b). ★ E-mail:sotillo@mpia.de Following a suggestion originally by Toomre (1977), galaxy mergers have also been shown to be able to drive morphological transformations in galaxies, i.e. to turn disky and rotationally-supported spiral galaxies into elliptical and dispersion-dominated ones. Already a few decades ago, it had been demonstrated via N-body simulations of isolated and idealized systems that the merger of two stellar disks generates descendants with spheroidal morphology (e.g. Barnes 1988;Hernquist & Barnes 1991;Barnes 1992;Hernquist 1993). Yet, more recently, it has also been found, with full-physics, cosmological hydrodynamical simulations of large samples of objects, that galaxies with similar stellar mass have undergone, on average, similar merger histories irrespective of their stellar morphology upon inspection (e.g. Rodriguez-Gomez et al. 2017;Martin et al. 2018;Jackson et al. 2020, with Illustris and Horizon-AGN).
this is the case not only in the low-redshift Universe (e.g. results with the SDSS and GAMA surveys, Park et al. 2007;Kelvin et al. 2014, respectively) but also up to ∼ 1.2 (e.g. results with VIMOS-VLT data by Ilbert et al. 2006).
Recent observations have enabled estimates of the merger histories of both the MW and M31. It is now possible to measure the trajectories of satellite galaxies surrounding the Galaxy and Andromeda, and, for the case of the MW, also the positions and velocities of individual stars (e.g. with RAVE, APOGEE, GAIA, Kunder et al. 2017;Majewski et al. 2017;Lindegren et al. 2016), both in the disk and the stellar halo, and in the satellite galaxies, together with ages and metallicities for many of these stars (e.g. with LAMOST and GALAH, Deng et al. 2012;De Silva et al. 2015). A series of stellar streams have been identified in the stellar halo that are associated with merger events.
Recently, Helmi et al. (2018) and Belokurov et al. (2018) have used data from the Gaia mission (Gaia Collaboration et al. 2018) and inferred that the Galaxy has undergone a major merger event sometime at = 1 − 2 (see also Bonaca et al. 2020;Xiang & Rix 2022): the merging galaxy has been dubbed Gaia-Sausage-Enceladus (GSE) and its mass at the time of the collision has been estimated to be ∼ 5 − 6 × 10 8 M in stars . The GSE merger has been associated theoretically with the formation of the inner stellar halo and the thick disk (Grand et al. 2020). Dynamical and chemical evidence of additional events in the MW's past has also been found, i.e. from the stellar remnants of past merging galaxies: these include Sequoia, accreted at a comparable epoch as GSE (Myeong et al. 2019); Thamnos (Koppelman et al. 2019); Helmi Streams (Helmi et al. 1999), that originated from a dwarf galaxy accreted approximately 5-8 Gyr ago (Koppelman et al. 2019), and the long-known Sagittarius Stream (Ivezić et al. 2000). The latter is still clearly identifiable also in configuration space and is being stripped from the Sagittarius dwarf of ∼ 5 × 10 8 M in stars, which is in the process of merging with the MW and whose infall time has been estimated to have occurred ∼8 Gyr ago (Dierickx & Loeb 2017). These, together with the ongoing mergers with the Large and Small Magellanic Clouds (with stellar masses of ∼ 2.7 × 10 9 M and ∼ 3.1 × 10 8 M , respectively (van der Marel 2006), and thus with a current stellar mass ratios of ∼1:10 and ∼1:100) constitute the six most massive accretion events across the known history of the MW, with the GSE being the largest in estimated mass ratio.
For M31, observational constraints are harder to obtain. However, for example, D'Souza & Bell 2018 recently suggested, by combining the results of cosmological models of galaxy formation and measurements of the ages and metallicities of halo stars, that Andromeda underwent a major merger about 2 Gyr ago (D'Souza & Bell 2018) and that the observed satellite galaxy M32 is the remnant core of the secondary galaxy involved in the merger: the latter produced the giant stellar stream that is located in Andromeda's stellar halo, generated a recent burst of star formation, but did not destroy the stellar disk.
The merger histories of MW-mass haloes have also been studied from a theoretical perspective, and via cosmological simulations. For example, Stewart et al. (2008) used DM-only or gravity-only N-body ΛCDM simulations and quantified that 95 per cent of MWmass haloes, i.e. gravitationally-collapsed objects of total mass of ∼ 10 12 M h −1 at = 0, have merged with at least one other halo more massive than 1:20 in the last 10 Gyr, this fraction reducing to 70 per cent for mergers with haloes more massive than 1:10. They concluded that halo-halo mergers involving (total) mass ratios of up to 1:5 must not destroy the stellar disks of the hosted galaxies if approximately two thirds of the observed galaxies are disky at = 0.
The reasoning mentioned above encapsulates and exemplifies a debate that has gone on for a few decades: namely, how to reconcile the fact that disk galaxies are found everywhere and at the same time the number of mergers that potentially would transform or destroy them are frequent. In the case of the Galaxy, Toth & Ostriker (1992) claimed via analytical arguments that a thin stellar disk similar to the Milky Way's is only compatible with no minor nor major mergers in the last 10 Gyr, a conclusion similar to that derived via observational data as early as by Wyse (2001) and Hammer et al. (2007). However, from a theoretical and numerical perspective, it has long been known (e.g. Hernquist 1989; Barnes & Hernquist 1996) that the formation and evolution of stellar disks -also and especially in the presence of perturbation mechanisms such as major mergerscannot be addressed without accounting for the role and availability of gas. Numerical dissipational simulations of idealized mergers by e.g. Naab et al. 2006;Robertson et al. 2006 showed that the remnants of gas-rich mergers can lead to the formation of a spiral and not an early-type galaxy.
In the full cosmological context, i.e. with zoom-in cosmological N-body and hydrodynamical simulations of galaxies, the formation of thin stellar disks in MW-mass galaxies could be achieved only after solving for the "angular momentum problem" or "overcooling catastrophe" and after introducing appropriate (stellar) feedback and star formation recipes (Guedes et al. 2011;Agertz et al. 2011) and refined numerical approaches Sĳacki et al. 2012;Vogelsberger et al. 2012). A number of successes followed, with simulations capable of returning large, disk-dominated galaxies that resemble the Milky Way in many respects (e.g. Martig et al. 2012;Stinson et al. 2013;Marinacci et al. 2014). Still, as of just some years ago, such zoom-in simulations of one or a few objects seemed to indicate that disk-dominated galaxies with stellar structures in accord with observations of the Galaxy could emerge, if, and only if, there is no major merger at 1 (Rix & Bovy 2013). The advent of large-volume cosmological hydrodynamical galaxy simulations ) such as Illustris (Vogelsberger et al. 2014a,b;Genel et al. 2014;Nelson et al. 2015) and EAGLE (Schaye et al. 2015;Crain et al. 2015) allows us to reconsider the formation of disk galaxies from major mergers in a quantitative manner, thanks to thousand-strong galaxy samples, spanning wide mass ranges and without prior constraints on the past assembly histories or environments. For example, with Horizon-AGN (Dubois et al. 2014) and Illustris galaxies, Martin et al. (2018) and Peschken et al. (2020) have confirmed that if enough gas is available during a merger event, new stars can form and can reform a stellar disk even if the latter gets destroyed in the collision. A similar scenario has been previously explored using zoom-in simulations by Sparre & Springel (2017), who considered four ∼MW-mass (at = 0) galaxies that were disky and underwent one major merger at ∼ 0.5. Focusing on very massive disk galaxies from Horizon-AGN ( * ≥ 10 11.4 M ), Jackson et al. 2020 showed that, in the majority of cases, pre-existing stellar disks either survive major mergers or reform shortly after a merger (a pathway in qualitative agreement with observations, e.g. Rothberg & Joseph 2004;McDermid et al. 2006) and, in fewer cases, stellar disks form after a spheroidal galaxy is established (Hau et al. 2008;Kannappan et al. 2009).
However, so far no quantitative analysis of the merger history and disk survival through mergers has been put forward with a focus on simulated analogues of the Galaxy and Andromeda. On the one hand, zoom-in simulation campaigns have begun to build up impressive samples of Milky Way-like galaxies: however, these are typically not unbiased in formation history/isolation (e.g. Auriga, Grand et al. 2017) or are selected in relatively narrow ranges of total halo mass (e.g. Artemis, Font et al. 2020). On the other hand, large-volume sim-ulation projects have lacked the resolution to capture stellar disks as thin as a few hundred parsecs. Large, unbiased samples at relatively high-resolution are now possible with the cosmological magnetohydrodynamical simulation TNG50 (Nelson et al. 2019b;Pillepich et al. 2019), the highest-resolution run of the IllustrisTNG project (www.tng-project.org). There we can identify approximately two hundred galaxies at = 0 with global properties similar to the MW's and M31's and with stellar disk heights as small as 100−200 parsecs. In addition to gravity and hydrodynamics, these simulations account for the effects of stellar and AGN feedback, magnetic field evolution, gas cooling and heating, etc.
In this paper, we hence describe the TNG50 cosmological simulation, the method to create the merger trees, and the selection criteria for the MW/M31 analogues in Section 2. We quantify the merger and assembly histories of TNG50 MW/M31-like galaxies in Section 3. In Section 4 we delve into the evolutionary pathways and physical mechanisms at play in the case of MW/M31 analogues that underwent at least one major merger in recent epochs, e.g. in the last 5 Gyr. We compare the = 0 structural and global properties of these galaxies with the rest of the TNG50 MW/M31-like objects in Section 5 and summarize and discuss our findings in Section 6.

The TNG50 simulation
In this paper, we focus on Milky Way and Andromeda-like (MW/M31-like) galaxies realized within the TNG50 cosmological simulation. TNG50 (Nelson et al. 2019b;Pillepich et al. 2019) is the highest-resolution run of the IllustrisTNG project (Naiman et al. 2018;Marinacci et al. 2018;Pillepich et al. 2018b;Nelson et al. 2018;Springel et al. 2018). It consists of a periodic cubic volume with side length of ≈ 51.7 comoving Mpc, contains 2160 3 dark matter (DM) particles and the same initial number of gas cells. The mass of the DM particles is uniform, DM = 4.5×10 5 M , and the average mass of the gas cells (and stellar particles) is baryon = 8.5 × 10 4 M . The code A (Springel 2010), which makes use of an unstructured moving mesh that is based on the Voronoi tessellation of the space and that is automatically adaptive in resolution, solves the gravity and magnetohydrodynamical equations in an expanding universe from ≈ 127 to = 0. Additionally, a series of physical processes are included to follow galaxy formation: these and their implementation are described in detail in the methods papers by Pillepich et al. 2018a;Weinberger et al. 2017. Crucially and differently from most high-resolution simulations of Milky Way-like galaxies (e.g. zooms without SMBH feedback by Wetzel et al. 2016;Buck et al. 2020;Renaud et al. 2021), the IllustrisTNG model not only includes star formation and its ensuing expected feedback and metal enrichment, but also magnetic fields and the seeding, growth, and feedback from supermassive black holes (SMBHs). All simulations in the Illus-trisTNG series adopt values for the cosmological parameters from Planck Collaboration et al. (2016).

Halo identification and histories of simulated galaxies
Throughout our analysis, the identification of haloes and subhaloes in the simulated domain is performed with the Friends-of-Friends (FoF; Davis et al. 1985) and S (Springel et al. 2001) algorithms, respectively. The latter identifies all particles and resolution elements that are measured to be gravitationally-bound to a self-gravitating structure.
We define as galaxies those subhaloes with non-vanishing stellar mass. The "central" galaxy of a FoF halo is typically the most massive one, whereas the remaining subhaloes within a FoF halo are called satellites. Throughout the analysis, we consider only subhaloes with cosmological origin, i.e. non-spurious subhaloes (Nelson et al. 2019b), and galaxies with at least 50 stellar particles.

Merger histories and merger summary statistics
Information for all resolution elements, subhaloes, and haloes is stored for one hundred snapshots, starting at = 20 (see Nelson et al. 2019b, for details). At recent epochs, the time spacing between consecutive snapshots is approximately 150 million years. The (sub)haloes at different snapshots are linked to describe their merger histories (or merger trees) using the S (Rodriguez-Gomez et al. 2015) and LH T  algorithms. Both return similar results. In this work, we use S G , which connects subhaloes across cosmic epochs based on their star particles and star-forming gas cells, instead of DM particles.
According to S G (described in detail Rodriguez-Gomez et al. 2015), the main progenitor of a galaxy is the one with the most massive history behind it: we will refer to its temporal evolution as the "main-progenitor history" (or "main-progenitor branch").
In the following, we describe the merger history of a galaxy by analyzing the time of its mergers and the stellar-mass ratio between the two involved galaxies, i.e. between the progenitors. Namely: • Time of a merger: As per Rodriguez-Gomez et al. (2015), the time of a merger ( merger ) is given by the earliest snapshot when both progenitors are identified by S as part of the same unique subhalo. Such a merger time represents the time of coalescence or just after it, i.e. at the next available snapshot.
• Time of the merger-ratio estimate: The stellar mass ratio (the smallest of the ratios between the secondary and the primary) is measured at the time when the secondary reached its maximum stellar mass ( max ). This can be considered as the time when a merger starts.
• Merger duration: The duration of a merger is the time elapsed between max and the time of coalescence, i.e. the time of a merger. This is therefore similar to the infall time, i.e. the time spent by the satellite in orbit within the host halo.
Mergers are classified into three types: major, whereby the stellarmass ratio is 1:4 or larger; minor, with ratios between 1:10 and 1:4, and mini (below 1:10). In the case of several mergers occurring at the same time, i.e. between the same pair of snapshots for the same galaxy, we account for all of them, separately, and their stellar mass ratios are recorded separately for each binary event.
In comparison to the default functioning of the S G code, we apply additional measures to avoid counting as mergers the interactions with spurious subhaloes or with subhaloes that are resolved with too few stellar particles, i.e. to account for the selections mentioned above. The former may be objects without a cosmological origin (e.g. fragments of disks or other types of galactic sub structures) that S identifies as subhaloes. They are identified with the SubhaloFlag (Nelson et al. 2019a) and are here excluded from the merger statistics. Galaxies whose main progenitor history is very brief (less than three snapshots before the time of coalescence) are also excluded. Finally, we also neglect all mergers between galaxies with fewer than 50 star particles (∼ 5 × 10 6 M ), across cosmic epochs. Thinking about the merger history of massive galaxies like MW/M31 analogues at = 0, this allows us to avoid counting as major those mergers whose max may occur very early on in the cosmic history of a galaxy, when primary and secondary objects might have had comparable masses, but whose actual coalescence occurs several billion years later and at times when the two progenitors actually differ by many orders of magnitude in mass.
The three conditions above allow us to obtain a clean and complete catalog of merger events, particularly at lower redshifts. On the other hand, these conditions also imply that our analysis is neither sensitive nor complete for the merger history of MW/M31-like main progenitors prior to ∼ 5.

In-situ vs. ex-situ stars
The cosmological simulations allow us to also follow the evolutionary tracks of individual stellar particles, so that we can identify the galaxy where each formed and the location at subsequent times. Based on the merger trees and for any galaxy at = 0, we categorize its stars as in situ if they formed in a progenitor that belongs to the main-progenitor branch of the galaxy, and ex situ otherwise, according to the definitions and catalogs described in Rodriguez-Gomez et al. (2016). Ex-situ stars are hence stellar particles stripped and consequently accreted into a galaxy from its mergers and orbiting satellites.

Galaxy and stellar-particle properties: circular orbits and diskyness or D/T ratio
Throughout the paper, the following conventions are intended unless otherwise stated.
Galaxy stellar mass: sum of the mass of all gravitationally-bound star particles.
Total host mass: the total mass of the FoF group in which a galaxy is found, 200c,Host , measured inside a sphere whose mean density is 200 times the critical density of the Universe (at the considered time).
Stellar half-mass radius of a galaxy ( 1/2 ): the radius of the sphere, centered at the particle with the minimum gravitational potential energy in the galaxy, that contains half of the stellar mass.
Star formation rate (total, SFR) is measured by adding the individual instantaneous SFRs of all gas cells that are gravitationally-bound to a galaxy at the time of inspection.
Other targeted galaxy properties are described in the course of the analysis.
For any given galaxy, the orbital properties of its stars are described via the circularity parameter, according to the definition by Scannapieco et al. (2009): = z / circ , where z is the specific angular momentum of the star in the direction perpendicular to the galactic disk, and circ is the specific angular momentum of a star at the same radius, but on a perfectly circular orbit. Namely, circ = circ , with circ =

√︁
(≤ )/ being the circular velocity of the galaxy at the considered radius. The up vector or vertical axis of a galaxy is chosen as the direction of the total angular momentum of all stars within two half-mass radii. Stellar orbits with 0.7 are considered circular, i.e. in rotational motion.
We define the diskyness or disk-to-total (D/T) ratio of a galaxy as the fractional mass of stars in circular orbits, i.e. with > 0.7, minus the fraction of stars with < −0.7, hence removing the contributions of bulge and stellar halo and assuming that the latter are symmetric around = 0. The stellar mass in circular orbits and the total stellar mass are evaluated, for the purposes of D/T, within an aperture of five stellar half-mass radii.
The exact choice of separating circular and cold orbits from the rest with the threshold > 0.7 is not crucial: Aumer et al. (2013) found that the D/T ratio obtained from our circularity definition is roughly equivalent to the values obtained with another usual definition of the circularities, based on the ratio of the specific angular momentum of the star and the maximum specific angular momentum at the specific binding energy of the star, using in both cases a similar circularity threshold.

Sample selection: MW/M31-like galaxies in TNG50
From the TNG50 box, which at = 0 samples hundreds of massive galaxies, we select the best analogues of the Milky Way and Andromeda based on their = 0 properties, i.e. without imposing conditions on their evolution or any restrictions on their morphological changes with time.
In practice, we use the selection presented and motivated in Aumer, White, Naab & Scannapieco (Pil) and already used by Engler et al. (2021); Pillepich et al. (2021) and others. To be selected as "MWor M31-like", a TNG50 galaxy must meet, at = 0, the following three conditions, based on galaxy stellar mass, stellar diskyness, and isolation: A) Galaxy stellar mass: log 10 (M * /M ) within [10.5, 11.2], with stellar mass measured within a 3D circular aperture of 30 kpc. B) Diskyness: either of the following conditions: (a) Stellar minor-to-major axis ratio (Pillepich et al. 2019), / , < 0.45, whereby and are respectively the minor and major axis of the ellipsoidal distribution of the stellar mass between one and two times the stellar half-mass radius. (b) Disky by visual inspection of 3-band images, both face-on and edge-on.
C) Isolation: no galaxy with galaxy stellar mass ≥ 10 10.5 M within 500 kpc distance and total host mass 200c,Host < 10 13 M .
There are 198 galaxies in the TNG50 simulation that match such conditions.
The stellar mass constraint is an envelope of the most accurate available estimates of the stellar mass of both the MW and M31 (e.g. Bland-Hawthorn & Gerhard 2016;Boardman et al. 2020). The diskyness criteria aim at including all simulated galaxies with a global stellar disk morphology and with spiral arms. Whereas typically, disky galaxies are defined as those with minor-to-major axis ratio < 0.33 (and middle-to-major axis ratio > 0.66; see e.g. Fig.  8 of Pillepich et al. (2019) and reference therein), the connection of this metric to other measures of "diskyness" (e.g. kinematic or photometric D/T ratios, scale length to scale height ratio, etc.) can be very varied, depending also on where, within a galaxy's body, the measurements are taken. Here we adopt a looser upper limit on the minor-to-major axis ratio to accommodate for the ratio of e.g. the stellar disk height to the stellar disk length of the Galaxy, particularly of its geometrically-thick component. The latter measures ∼ 0.45 (in comparison to ∼ 0.11 for the geometrically-thin disk) when adopting the best values for the structural properties of the Milky Way disk from Bland-Hawthorn & Gerhard (2016, thin and thick disk lengths = 2.6 ± 0.5 and 2 ± 0.2 kpc, respectively, and thin and thick   Figure 1. Mass assembly histories of MW/M31 analogues in TNG50, for galaxy stellar mass (left) and total host mass (right). We show the complete TNG50 MW/M31 sample in black and a sub-sample of galaxies with recent major merger in blue (see next Sections): median tracks at fixed cosmic time are marked as thick curves, whereas shaded areas denote interquantiles (10th-90th percentiles). Additional curves represent results from other simulations (colored lines) or observationally-derived constraints (grey lines with empty markers). For TNG50 and the additional simulations, the subhalo's total stellar mass is plotted. For = 0 MW/M31-like galaxies, 10 per cent of their stellar mass is assembled by = 2 and 50 per cent by = 1. At = 2 the scatter in galaxy stellar mass extends for two orders of magnitude (1.8 × 10 8 − 2.2 × 10 10 M , 10th-90th percentiles), in contrast to a = 0 selection ranging 0.7 dex. In terms of total host mass, = 0 MW/M31-like galaxies have assembled 50 per cent of their halo mass by = 2. We show annotations with observational estimates for the stellar mass at = 0 of the MW (in magenta, as an envelope of the values from Flynn et al. 2006;Licquia & Newman 2015;Bland-Hawthorn & Gerhard 2016) and M31 (orange, from Geehan et al. 2006;Barmby et al. 2006;Tamm et al. 2012;Sick et al. 2014). disk heights at the Sun location = 300 ± 50 and 900 ± 180 pc). The isolation criterion avoids the presence of a galaxy, with mass equal to the lower MW estimate or larger, at a distance closer than 500 kpc -the distance between the Galaxy and Andromeda is ∼ 770 kpc (McConnachie et al. 2005;Riess et al. 2012). Additionally, the requirement on total host mass is meant to exclude galaxies located within a very massive galaxy group or cluster, as we know this is not the environment of the Galaxy and Andromeda. However, it is more relaxed than requiring a galaxy to be the central of its halo, while allowing certain galactic environments such as those similar to the Local Group, whereby the Galaxy and Andromeda may be sharing the same dark matter halo.
As the galaxies are selected from a volume-limited sample, a larger number of TNG50 objects have stellar mass closer to the MW's than to the more massive M31's. In fact, among the 198 TNG50 MW/M31 analogues, 130 galaxies have masses more compatible with the MW's (≤ 10 10.9 M ), and will be dubbed 'MW-mass' and 68 are instead more representative of Andromeda and will be dubbed 'M31-mass'.

Mass assembly
We quantify the assembly history of TNG50 MW/M31-like galaxies by plotting their mass across cosmic epochs, i.e. along their individual main progenitor branches. The evolution of the galaxy stellar mass (left) and of the total host mass (within 200c , right) of all MW/M31 analogues are shown in Fig. 1. For TNG50 and the additional simulations the total (i.e. gravitationally-bound) stellar mass is plotted. Thick curves denote the median evolutionary tracks across the galaxy populations at fixed redshifts, whereas the shaded areas indicate the 10th-90th percentiles. It should be noted that, although the medians appear smooth, individual assembly histories may exhibit assembly histories that are very different from the median curves. The selection by galaxy stellar mass at = 0 -in the range of 3.2×10 10 −1.6×10 11 M , within 30 kpc aperture as per Section 2.4, left panel -translates into a = 0 host mass range of M 200c = 8.3 × 10 11 − 2.5 × 10 12 M at = 0, i.e. 0.48 dex, within 10th-90th percentiles (right panel): this is a consequence of the effective stellarto-halo mass relation resulting in TNG50 (Engler et al. 2021). On average, the host haloes of MW/M31-like galaxies grow by a factor of ∼ 7 between = 3 and today, while the average galaxy stellar mass of their progenitors at = 3 was a factor of ∼ 40 lower than today.
The progenitors of MW/M31-like galaxies span a host mass range (within 10th-90th percentiles) at = 2 similar to that at = 0, namely 0.5 − 0.6 dex: on the other hand, the progenitors of MW/M31-like galaxies at e.g. = 2 span almost three orders of magnitude in stellar mass across all sampled galaxies: 1.8 × 10 8 − 2.2 × 10 10 M In the tables, different rows correspond to different stellar mass ratios. Top panel: fraction of galaxies that have undergone at least one merger in different time periods (lookback time measured from = 0) and different stellar mass ratios: major (red), minor (blue), minor or major (purple), mini (green) and any (black). Dots represent the exact periods for which the measurements are available. Bottom panels: fraction of galaxies undergoing different types of mergers for three defined time periods. Left: last 5 Gyr ( ∼ 0.5). Middle: since = 1 (lookback time ∼ 7.98 Gyr). Right: since = 2 (lookback time ∼ 10.51 Gyr). The color denotes the fraction of galaxies undergoing a certain number of mergers. These fractions are not cumulative. 18 per cent of the MW/M31-like galaxies have recent (i.e. over the last 5 Gyr) major mergers. 34 per cent have a major merger since = 1.
These theoretical predictions are of the essence to connect galaxy populations across cosmic time for the purposes of contrasting e.g. MW-like galaxies with their expected progenitors at high redshifts. For comparison, in Fig. 1 we show the corresponding results from previous cosmological simulations of MW-like galaxies, albeit differently selected: solid and dashed thin lines. The thirty zoom-in Auriga MW analogues (yellow dashed curve, for their median growth, Grand et al. 2017), where galaxies are selected as isolated haloes at = 0 with 200 in the range (1 − 2) × 10 12 M , have a very similar median stellar assembly history to those from TNG50. For the Illustris simulation (Torrey et al. 2015a), where MW analogues are all galaxies with stellar mass within the range (4 − 5) × 10 10 M , the median stellar mass is lower at all redshifts (red solid curve). The six NIHAO-UHD MW-like galaxies (green dashed curves, Buck et al. 2020) span a similarly wide range in mass growth across all cosmic epochs as TNG50. On the other hand, the Eris galaxy (orange dashed curve, Guedes et al. 2011) exhibits a much earlier mass growth (Pillepich et al. 2015), similarly as to the FIRE-2 MW-mass zoom-ins by Garrison-Kimmel et al. (2018), of whom we show the growth curves for the least and most massive at = 0 (dotted cyan).
In Fig. 1, assembly histories obtained by assuming that galaxies preserve their number density in time (e.g. van Dokkum et al. 2013) are shown for contrast. As in the original Illustris simulation (Torrey et al. 2015a,b), also according to TNG50 and for MW/M31-like galaxies only, we find that the average stellar mass evolution inferred via a constant comoving number density assumption is systematically shallower than when tracking galaxies via their merger trees, with ∼ 3 progenitors characterized by galaxy stellar masses a factor of a few larger than those obtained along the main progenitor branches of galaxies. Therefore, the mass evolution from constant comoving number density assumptions cannot be used to validate the results of zoom-in simulations (Buck et al. 2020).  . Stellar-light composite images of the 31 MW/M31-like galaxies from TNG50 at = 0 that have undergone a recent major merger (i.e. in the last 5 Gyr). Every disk galaxy is shown in face-on and edge-on projections. Each panel spans horizontally 80 kpc, and each galaxy is identified by its stellar mass and Subhalo ID number in the TNG50 = 0 catalog. These galaxies exhibit a great diversity in stellar sizes and in terms of stellar structures within the disks (see text for comments on galaxies with Subhalo IDs 400973 and 41961): for example, there are cases with marked spiral arms and others with central bars. Fig. 2 shows the merger statistics of MW/M31-like galaxies according to TNG50, selected as described in Section 2.4. To our knowledge, this is the first time that such statistics, based on galaxy stellar masses rather than halo masses, are quantified and are obtained from a cosmological hydrodynamical simulation where MW/M31-like galaxies are selected based on observable (rather than host halo) properties.

Merger statistics
The upper panel shows what fraction of MW/M31 analogues undergo at least one merger over varying past periods of time (in lookback Gyr from = 0), for different stellar mass ratios: major, minor, major or minor, mini and any ratio. In the lower panels, we give the fractions of galaxies that experienced varying number of mergers, over three time periods: from left to right, in the last 5 Gyr (i.e since z ∼ 0.5), since = 1 (lookback time ∼ 7.98 Gyr), and since = 2 (lookback time ∼ 10.51 Gyr). The color code in the lower panels represents the fraction of MW/M31-like galaxies that underwent a certain number of mergers -these fractions are not cumulative. In all panels, mergers are counted based on the times of coalescence, i.e. the time of the mergers as defined in Section 2.2. About thirty per cent of MW/M31-like galaxies have undergone at least one major or minor merger (i.e. at least one merger with stellar mass ratio larger than 0.1) over the last 5 billion years: this fraction increases to 48 (72) per cent since = 1 ( = 2).
Interestingly, for MW/M31-like galaxies, on average, the frequency of past major mergers is similar, if not slightly larger, than the frequency of minor mergers. Not only are the fractions of galaxies undergoing at least one minor or one major merger similar (upper panels), but also comparable are the fractions of galaxies undergoing different numbers of minor and major mergers since a fixed period of time (lower panels).
Whereas it has often been assumed, also prior to the newest results with Gaia, that our Galaxy has had a very quiet merger history at least since ∼ 1 with no major mergers since then (see Introduction), according to TNG50 and our selection, the fraction of MW/M31-like galaxies that have merged with at least one other similarly-massive galaxy since ∼ 1 is about 30 per cent. This fraction is slightly lower, 27 per cent (35 of 130 galaxies), if we only consider TNG50 analogues with stellar mass below 10 10.9 M , i.e. if we exclude Andromeda-mass galaxies.
From the lower panels of Fig. 2, it can also be seen that 9 per cent of MW/M31-like galaxies have not merged with any other galaxy, irrespective of its stellar mass, over the last 5 billion years. On the other hand, a handful (12 per cent) of the selected TNG50 galaxies have experienced multiple major mergers since = 2, as many as between 3 and 5. Yet, more than one merger event per galaxy (major or minor) since e.g. = 1 is infrequent 1 .
It is important to notice that the numbers summarized in Fig. 2 are based on a full-physics (i.e. not DM only) model for galaxy formation in the full cosmological context, they apply to a specific galaxy selection, and the merger mass ratios are characterized in terms of the stellar mass (and not DM mass) of the merging objects. Fig. 2 can be considered as the stellar-mass based update of e.g. Figs. 5 and 6 of Stewart et al. 2008, which were based on N-body only simulations, and the TNG updates of the merger estimates for MW-mass haloes given for Illustris by Rodriguez-Gomez et al. 2015 and quoted in the Introduction.

MERGERS AND DISK SURVIVAL
The findings in the previous Section imply that a selection of TNG50 galaxies at = 0 have global properties -stellar mass, stellar diskyness, and large-scale environment -similar to those of the Galaxy and Andromeda and yet have undergone at least one major merger as recently as over the last 5 Gyr. In particular, 31 of the 198 MW/M31 analogues of TNG50 (16 per cent) are found in this subsample and are hence the focus of the rest of the paper. The choice to focus on galaxies with the time of their last major merger at 5 billion years ago, instead of e.g. 6 or 4 billion years ago, is somewhat arbitrary but is simply meant to qualitatively encapsulate the phenomenology of relatively recent major mergers. We will drop the binary classification whenever instructive and possible.

MW/M31-like galaxies as survivors of recent major mergers
In Fig. 3, we show stellar-light composite images at = 0 for the subsample of TNG50 MW/M31 analogues with a major merger within the last 5 billion years, in face-on and edge-on projections. But for a couple of cases -Subhalo IDs 400973 and 419618, for which the identification as MW/M31 analog is in fact borderline -they pass the morphological selection criterion based on the shape of the inner stellar mass distribution but have D/T mass ratio of ∼ 0.2 -, these galaxies exhibit clear stellar disk morphologies, although with a great variety in structural properties and extents. Their median D/T mass ratio at = 0 (as defined in Section 2.3) is 0.40, but ranges between ∼ 0.12 and ∼ 0.78 across the whole sample of Fig. 3. For comparison, the D/T ratios of the full MW/M31-like sample range from ∼0.10 -0.90 with a median of 0.55 and six systems below 0.1, all with counter-rotating structures (see also Joshi et al. 2020).
The geometrically-thin stellar disks of MW/M31-like galaxies with recent major mergers have sech 2 heights measured at 4 times the disk scale length spanning between ∼50 pc and 2.9 kpc, with median and average of 0.93 and 1.14 kpc (see Section 5 for more details). In numerous cases, the stars are organized in very prominent spiral arms and grand-design spiral systems, which are sometimes distributed asymmetrically. There are also central bars with an extension of ∼ 1 − 4 kpc (e.g. Subhalo IDs 522983, 523548 and 540920). From the edge-on views, we can appreciate that the stellar haloes are also diverse, even among galaxies with similar stellar disk size: from faint haloes of stars that barely extend beyond the disk, to others that are appreciable in the images for up to a few tens of kpc above and below the disks. A quantification of the mass in the different stellar components (disks, bulges, and haloes) will be given in Section 5.
In Fig. 4, we provide additional statistics and properties of the last major mergers experienced by the MW/M31-like galaxies with a recent major merger (blue), in comparison to those of the whole MW/M31-like sample. Firstly, as appreciable also from Fig. 2, 15 per cent, i.e. 30 galaxies among the TNG50 MW/M31 analogues have never undergone a major merger, i.e. never since ∼ 5, see Method Section for details. In the top left panel of Fig. 4, these are reported as an orange bar: 6 of those are M31-mass galaxies.
If we look at the times when the last major mergers of MW/M31like galaxies occurred (time of coalescence, top left panel of Fig. 4), we see that 83 galaxies (42 per cent) experienced their last major merger more than 10 Gyr ago, whereas 95 galaxies (48 per cent) had their last major merger more than 9 Gyr ago -this compares to ∼35 per cent for * ∼ 10 10 M disk galaxies according to Font et al. (2017)). Importantly, the MW/M31-like galaxies with a recent major merger are not particularly biased in their merger times and span the entire final 5 billion years of cosmic evolution. Magenta and orange shaded vertical bands denote the current estimates of the last major mergers of the Galaxy, namely of GSE (Helmi et al. 2018;Myeong et al. 2018;Chaplin et al. 2020), and of Andromeda (D'Souza & Bell 2018), respectively.
In the top right panel, we see that the durations of the last major mergers (see Section 2.2) span from ∼ 0.04 to ∼ 3.5 Gyr. The longest mergers are more frequently recent, owing to the shorter dynamical times when the Universe was younger: 77 per cent of the mergers lasting longer than 1.5 Gyr occurred within the last 5 billion years. On average, TNG50 MW/M31-like galaxies experienced major mergers that carried on for 0.7 Gyr (median of 0.5 Gyr), this duration increasing to 1.4 Gyr for those happening since 0.5 (median of 1.3 Gyr). Dividing the TNG50 MW/M31-like sample according to stellar mass, and considering only the galaxies that underwent at least one major merger over their history, the last major merger of MW-mass (M31-mass) galaxies occurred on average about 10 (9) billion years ago. The median merger duration is approximately 0.61 and 0.46 Gyr for the two sub-samples.
The stellar mass ratios of the merging galaxies are quantified in the bottom panels of Fig. 4, evaluated at max (the time of the maximum stellar mass of the secondary, see Section 2.2). All ratios in the bottom left panel are above 0.25, by our definition of a major merger. The mass ratios closer to 1 are more frequent for older mergers (see our definitions of merger mass ratios in Section 2.2). Only one of the 31 recent major mergers has a mass ratio larger than 0.75. Consistently with the hierarchical model of galaxy assembly, the most recent major mergers have in general lower stellar-mass ratios: galaxies at recent times are more massive but massive galaxies that could merge with them are rarer. The absolute values of the galaxy stellar mass of the merging progenitors are shown in the bottom right panel of Fig. 4. MW/M31-like galaxies with recent major mergers (blue symbols) occupy the highest part of the range, at fixed stellar mass, and are compared to the progenitors of all MW/M31-like galaxies with the last major merger happening at any time (light blue symbols): the former subset of galaxies increases their mass by a small factor in the last few Gyrs, so we do not find galaxies with a recent major merger where both progenitors have low stellar mass and the descendant . Characteristics of the last major mergers of MW/M31-like galaxies in TNG50 at = 0: major mergers. In all panels, MW/M31-like galaxies that experienced their final major merger in the last 5 Gyr are depicted in blue; MW/M31 analogues with their last major merger at any cosmic time are indicated in light blue. Top left: Time of the last major merger (in lookback Gyr) for each of the 198 galaxies, in bins of 1 Gyr: 30 MW/M31-like galaxies did not undergo any major mergers (orange bin); 95 galaxies (approximately 50 per cent) underwent their last major merger more than 9 Gyr ago ( 1.3); 31 galaxies have had a late major merger, i.e. as recent as in the last 5 Gyr. The occurrence of the latter is distributed in an approximately uniform way in the considered period of time. Top right: Duration of the major mergers, defined as the time elapsed since the secondary reached its maximum in stellar mass ( max ) and the moment of coalescence -bins span 330 million years. 77 per cent of the major mergers lasting longer than 1.5 Gyr are recent. Bottom left: Stellar-mass ratios between the secondary and the primary, at max . Bottom right: Stellar mass of both progenitors, at max , versus the stellar mass of the galaxy at = 0. For the Galaxy and Andromeda, we show observational estimates with magenta and orange annotations, respectively: the mass-ratio and the time of the known last major merger are taken for the case of GES (Helmi et al. 2018;Gallart et al. 2019;Naidu et al. 2021) and for M32 (D'Souza & Bell 2018). grows then rapidly enough to be included in our MW/M31-like mass cut. Magenta and orange vertical shaded areas denote the current stellar mass constraints of the Galaxy and Andromeda, respectively: according to TNG50, even galaxies with lower, MW-like mass may experience a recent major merger and still be disky at = 0.
The figures above show that, according to TNG50, MW/M31-like galaxies with recent major mergers have interacted with relatively massive companions for substantial amounts of times, i.e. on average for ∼ 1.4 Gyr with secondaries of * ≈ 2 × 10 10 M , resulting in mergers of median stellar mass ratio of 0.41 (i.e. 1:2.5).

Gas availability during the mergers
But how is it possible that, despite having undergone major mergers as recently as in the last 5 billion years, galaxies can nevertheless exhibit marked disky stellar morphologies at = 0?
Previous studies (Hopkins et al. 2009;Stewart et al. 2009;Hoffman et al. 2010) had shown with idealized simulations of mergers that the gas content in the merging progenitors is a determinant factor for the outcome of a merger -see Introduction: it both conditions how destructive a merger is for the stellar component of the involved galaxies and how actively stars can form in the descendant galaxy. In Fig. 5, we hence analyze the gas availability during merger events for all the TNG50 MW/M31 analogues undergoing major mergers across cosmic epochs. Firstly, we estimate the gas content of each galaxy by accounting for all the gravitationally-bound gas at the time of the merger: this choice is supported e.g. by Sparre et al. (2021), who show with simulations that the gas in the circumgalactic medium and even in the outer halo can contribute to the star formation after major merger events between = 0.3 − 0.8 that produce disky, MWlike galaxies. Secondly, we compare the gas mass fractions of each system at coalescence (i.e. the gas mass over the total stellar mass) with that of central galaxies that, at the corresponding major-merger time of TNG50 MW/M31-like galaxies, have the same stellar mass (within a range of ±0.1 dex) as the galaxy that has resulted from the two progenitors. These galaxies that serve as reference are plotted in gray, with the median gas fractions marked with a gray solid line and the narrower and the broader shaded regions representing, respectively, the 25th-75th and the 5th-95th percentiles. MW/M31 analogues are marked as blue symbols and a median is added to facilitate comparison.
Two facts are noticeable. First, as anticipated (see also e.g. Pillepich et al. 2019 for TNG50 results) and as it is expected from observations of high-redshift galaxies, the average mass fraction of gas in galaxies decreases with cosmic time -for the control sample, the gas mass fraction drops by ∼1 dex from ∼ 3 (about 2 Gyr after the Big Bang) to ∼ 0.5 (about 5 Gyr ago). Generally, the MW/M31 analogues show a gas fraction, on average, comparable to that of the control sample of central galaxies in the similar mass range. Second, and importantly, the difference between samples increases slowly . Gas mass fraction ( gas / stars ) of TNG50 MW/M31-like galaxies at the time of their last major merger (time of coalescence). For comparison, in gray, we show the gas fractions of central galaxies in the same stellar mass range (±0.1 dex) at the corresponding merger times: the median is the solid curve. MW/M31-like galaxies with recent major mergers (blue symbols) have lower gas-to-star mass fractions, in absolute terms, than those with more ancient mergers -i.e. the gas mass fraction of galaxies decreases with time. However, they are more gas rich (0.1 − 0.2 dex) than the average central galaxies in TNG50. This gas richness can also be appreciated in the inset histogram, where we show the logarithmic difference of the gas mass fraction, for each of the galaxies, with respect to the median of the control sample in the considered time of the last major merger.
with cosmic time: we think that this is the very manifestation of the fact that galaxies that are selected to be disky at = 0 while having undergone a recent major merger constitute somewhat a biased subset. Namely, MW/M31-like galaxies with recent major mergers exhibit somewhat higher gas fractions (∼ 0.1 − 0.2 dex) than the average of other non-merging galaxies with similar mass at the corresponding epoch: 23 of 31 of them lie above the median value of the control sample. Whereas the gas-mass difference is not large, it appears sufficient to trigger star formation in TNG50 galaxies during the merger events, as we show next.

Star formation bursts triggered by gas-rich major mergers
The availability of gas enables star formation during the recent major mergers of TNG50 MW/M31-like galaxies.
We have examined the SFR evolution of the TNG50 MW/M31 analogues with recent major mergers and find that their SFRs can increase substantially in correspondence to their major merger events. Such SF bursts may occur at the time of coalescence or shortly after, and sustained SF may be in place in the descendants through = 0. In fact, SF bursts may happen even before the merging galaxies coalesce, specifically at the close pericentric passages of the secondary progenitor in its approaching orbit. About 24 of the 31 MW/M31 analogues with recent major mergers show appreciable bursts of SF triggered by the merger, seven of which are the depicted in Fig. 6.
There we show the evolution of the instantenous SFR along the main progenitor branch of selected galaxies (magenta) as a function of cosmic time (starting at 8 Gyr after the Big Bang), compared with their distance (in physical kpc) to the secondary progenitor involved in their last major merger: green curves. For each galaxy, i.e. row, the time of merger is represented with a vertical black solid line, whereas the last pericentric passages are distinguishable as the minima in the . Bursts of star formation in a selection of MW/M31-like galaxies in TNG50 with recent major mergers. SFR in the main progenitor (magenta) and distance between both progenitors (green) are plotted for seven example galaxies, one panel each. The time of the last major merger is marked with a solid vertical line; the last pericentric passages of the secondary progenitor are marked with dotted vertical lines. The galaxy mergers and the close galaxygalaxy interactions prior to coalescence can trigger substantial bursts of star formation. The SFR after the merger varies depending on the galaxy. distance curves: dotted black vertical lines. It can be clearly seen that, not only coalescence, but also close galaxy-galaxy interactions prior to the merger time can trigger important events of star formation in the main galaxy and consequently can plausibly alter its structure, morphology, and stellar mass content. It was previously shown that highly "bursty" SF was suppressed in galaxy major mergers of the original Illustris simulation due to its insufficient numerical resolution: in particular, this was demonstrated via zoom-in simulations of selected merger events with 40 times better mass resolution (Sparre & Springel 2016). Fig. 6 shows that the numerical resolution of TNG50 is sufficient to capture the compression of gas possibly due to the galaxy-galaxy interactions, the funnelling of gas towards the galaxy centers, higher gas density and hence to reproduce bursts of star formation triggered by mergers and galaxy interactions.

The cases of disks destroyed during the mergers and reformed vs. those surviving during the merger
By inspecting the time evolution of the SFR of the main progenitors of each MW/M31 analogue with recent major mergers and the time evolution and distributions of the orbital properties of their stars (as in Figs. 8 and 9), we find that two main scenarios or pathways are followed by the 31 TNG50 galaxies to be disky at = 0 after undergoing a recent major merger. These are as follows: (i) in 18 cases (58 per cent), the galaxy's stellar disk, being in place prior to the time of the merger or prior to the pericentric passages, is destroyed by the merger; nevertheless, the galaxy by = 0 reforms a stellar disk; (ii) in 11 cases (35 per cent), the stellar disk is not completely destroyed by the merger and therefore the galaxy retains its disky morphology down to = 0; Two remaining galaxies of the sample are difficult to classify. In one case (TNG50 Subhalo ID 419618 at = 0), the D/T mass ratio is actually quite low both prior to and after the merger as well as at = 0: ∼ 0.1 − 0.2 -we noticed this galaxy in Fig. 3 and possibly this should not have been included in the TNG50 MW/M31-like sample in the first place. In the case of TNG50 Subhalo ID 400973 at = 0, the D/T mass ratio around the time of the merger is biased low by the fact that (possibly another) merger event has triggered the production of new stars in counter-rotating orbits, making the study of its evolution rather specific. In both cases, the galaxies have experienced a major merger as recently as in the last one or two snaphots, i.e. in the last 150 − 300 million years, so much so that stellar shells are somewhat visible in their stellar light maps.
We have also explicitly checked whether cases exist whereby the recent major merger occurs at a time when the progenitor of the selected MW/M31-like galaxies does not exhibit yet a clear disky stellar morphology (or not anymore due to a previous merger). But for the two galaxies described above, among the 31 TNG50 MW/M31 analogues, this does not occur.
The frequency of major merger events in MW/M31 analogues and the different evolutionary pathways that the TNG50 simulation has allowed us to uncover are summarized in Fig. 7. In the various boxes from top to bottom we show the number of TNG50 MW/M31-like galaxies undergoing certain conditions and paths.
Somewhat similar scenarios had been found in previous studies, albeit without a specific focus, and hence selection, on the MW/M31 mass scale. For example, Jackson et al. (2020) used the Horizon-AGN cosmological simulation and reported about two main pathways for the existence of massive galaxies with disky stellar morphologies at = 0 ( * ≥ 10 11.4 M , above our upper end): surviving disks (30 per cent), owing to anomalously quiet merger histories, and spheroidal galaxies that form a young stellar disk (70 per cent) as a consequence of star formation triggered by a recent gas-rich mergers. According to Font et al. (2017) and the GIMIC simulation, on the other hand, stellar disks reform for one third of the * ∼ 10 10 M galaxies undergoing a major mergers in the last 7 Gyr. Also, Sparre & Springel (2017) showed that four * ∼ 10 10.5−11 M disk galaxies that were disky before their last major merger (at ∼ 0.5) were able to reform the disk after the merger destroyed it.
Figs. 8 and 9 exemplify the time evolution of one prototypical example galaxy for each of the two scenarios identified for TNG50 MW/M31-like galaxies. Plotted as a function of cosmic time, from top to bottom, are (more details in the caption): stellar and gas mass for the primary and secondary progenitors, SFR of the main progenitor, distance between progenitors, diskyness (D/T ratio) and the Figure 7. Summary of the evolutionary pathways uncovered via the TNG50 simulation for a sample of galaxies selected at = 0 as MW/M31 analogues. MW/M31-like galaxies may have had at least one major merger throughout their history or not at all: 85 and 15 per cent, respectively -second panels from the top). Among the former, third layer from the top, 31 galaxies, i.e. 18 per cent of the TNG50 MW/M31-like sample with at least one major merger ever, actually underwent one or more major mergers as recently as in the last 5 billion years. For most of them (58 per cent), the previously-existing stellar disk was destroyed but reformed; for another substantial fraction (35 per cent), the previously-existing stellar disk survived through the merger event. Inset diagrams show the evolution of D/T with time, centered at the time of the merger. Thin (thick) curves represent individual galaxies (medians). For the destroyed disk the drop of diskyness is much more noticeable. Two of 198 galaxies cannot be placed in either scenario -see text for details. In each box, we report the number of galaxies in absolute and relative values, being the latter calculated within the downstream sub-samples.
distributions of the orbit circularities of the stars formed in situ. Finally, images of the stellar mass density of the main progenitor are shown from face-on and edge-on projections, at the three snapshots immediately prior to the last major merger, at the snapshot after the merger and at = 0 -which helped in the selection of the galaxies in the MW/M31 sample. In each figure, the vertical solid, thick line represents the time of the last major merger, i.e. the time of coales-  . Evolutionary tracks of an example TNG50 MW/M31-like galaxy whose last major merger does not destroy the stellar disk. Panels as in Fig. 8. cence. The blue thin vertical lines in the top panels denote the time when the secondary reaches its maximum stellar mass and the stellar mass ratio is evaluated.
In the first pathway, i.e. Fig. 8, a disk of young stars grows in the period of a few billion years between the last major merger and = 0.
In these cases, the amount of gas available to form new stars is a key factor. From the Figures, it is also clear that, as already alluded to in Fig. 6, not only does the major merger affect and alter the galaxies, but also the close pericentric passages of the secondary progenitor can have a large impact on the final outcome of the merger and on . Connection between ongoing star formation and galaxy stellar morphology, i.e. diskyness, for TNG50 MW/M31-like galaxies with recent major mergers. We show the evolution along the main-progenitor branch of individual galaxies of their SFR (fuchsia) and of their D/T ratio (blue), as a function of cosmic time. We show random examples of galaxies whose stellar disk had been destroyed by their last major merger and regenerated (left column: 6 examples of 18 cases) and of galaxies whose stellar disk has not been destroyed (right column: 6 examples of 11 cases). The time of the last major merger is marked with a black solid vertical line; the times of the last pericentric passages are marked with black vertical dashed lines. The destruction of the stellar disks in the left column is identifiable in those cases where a galaxy undergoes a sudden drop in diskyness, and usually this occurs in the period between the last pericentric passage and the time of the merger. Compared to the destroyed disks, the galaxies on the right column exhibit milder and gentler drops in diskyness. In all cases, sustained star formation is in place during and after the merger, often with bursts at pericentric passages (see also Fig. 6) and is temporally coincident with a steady D/T increase after the merger and towards = 0, particularly so in the galaxies where the merger destroyed the disk. star formation. In fact, but for periods around coalescence (Fig. 8) and possibly at pericentric passages, new stars are born in circular orbits, i.e. with circularity at birth typically close to unity -this is the case at all times and for all galaxies (Pillepich et al. 2019). On the other hand, the current circularity of the in-situ stars that are in = 0 MW/M31-like galaxies may be very different, i.e. hotter, than that at birth, particularly for stars that formed at early times: top vs. third circularity panels.
For the galaxy of Fig. 8, which is an example of a last major merger that destroys the stellar disk of the main progenitor, the merger produces a drop in the stellar circularities of the galaxy's stars (top and bottom circularity panels). However, stars formed after the merger are born mostly in circular orbits, so that a new stellar disk is present at = 0. This sequence is also visible in the stellar density images. Fig. 9 represents a case of a disky main progenitor whose stellar morphology is not affected by the major merger: compared to the previous case, the secondary progenitor approaches the main progenitor more progressively, with multiple pericentric passages prior to the final merger.

Connection between in-situ star formation and diskyness during and after the mergers, and on the accreted mass
As pointed out above, within the TNG50 model, the orbits of stars at the time of formation are generally circular -at least at 2 − 3, see 7th panels from the top of Figs. 8 and 9. Namely, if star formation occurs, stars naturally form in circular orbits, because they originate from gas that is in rotationally-supported and disky configurations (see also findings by Pillepich et al. 2019). Here we further expand on the connection between star formation and galaxy stellar morphology.
In Fig. 10, we examine several examples of TNG50 MW/M31like galaxies with recent major mergers and for each we show the evolution in time of the instantaneous star-formation rate (fuchsia) and of diskyness D/T (blue) along the main progenitor branch. In these plots the time of the merger is marked with a vertical solid black line and the moments of the last pericentric passages with vertical dashed black lines. We show 6 random examples among the cases whose last major merger destroyed the disk but the latter reformed (left panels) and 6 random examples among those galaxies whose last major merger did not destroy the disk. stars formed BEFORE last major merger stars formed DURING last major merger stars formed AFTER last major merger . Connection between star formation, orbital circularities of the stars and the amount of in-situ star formation triggered by the merger and of the ex-situ stars brought in by the mergers. Top: circularities of the stellar orbits at the time of formation vs. at = 0, stacked for all 31 TNG50 MW/M31-like galaxies with a recent major merger, for stars formed before (before max ), during (between max and merger ), and after the last major merger, from left to right. In all panels, stars in perfectly circular orbits have a circularity = 1 (or −1 if they are counter-rotating); stars in radial orbits have a circularity of ∼ 0. Histograms show the total stellar mass per galaxy born in the respective periods of time. Bottom left: Fraction of in-situ stellar mass formed in the last 5 Gyr as a function of total stellar mass at = 0. Bottom right: Ex-situ stellar mass brought during the last major merger vs. in situ stellar mass in the last 5 Gyr. In both cases, TNG50 MW/M31-like galaxies with a recent major merger (blue) are contrasted with those with at least one major merger across their history but irrespective of when (light blue). Most of the stars were born in circular orbits, this is generally the case both before and after the times of the last major merger and, to a lesser extent, during the last major merger; however, the orbits of the stars born after the last recent major merger remain mostly unaltered, i.e. in circular orbits, all the way to = 0. Furthermore, recent major mergers also seem to trigger more in-situ star formation in the resulting galaxy than in galaxies with more ancient mergers, in addition to bringing large amounts of ex-situ stellar mass.
First, in the left panels, a noticeable feature is the decrease of the diskyness of the galaxy around the time of the merger or some million years before the merger: in the latter cases, this is due to the effects of the proximity of the secondary progenitor at the pericentric passages. The drop of diskyness is roughly at least one third of the D/T value before any interaction with the secondary. This is not the case for the galaxies in the right panel, where the changes in circularity fraction, D/T, are less pronounced.
Second, the evolutionary tracks of Fig. 10 confirm that, after the merger, a disk of young stars can form again: this can be appreciated in the fact that the D/T ratio increases progressively without new drops (the galaxies can undergo additional minor or mini mergers, but no major ones). A progressive increase in the diskyness is here correlated with sustained phases of star formation.
The connection between ongoing and recent star formation and stellar morphology is quantified and extended to the entire galaxy sample in Fig. 11. In the top panels, we compare the circularities of the stellar orbits at birth and at = 0, for the stars formed before, during, and after the last major merger: namely, before max , between max and merger , and after merger , i.e. after coalescence, respectively. Results are shown by stacking the orbital properties of all the stars in all the TNG50 MW/M31 analogues with a recent major merger -therefore the panel on the top left (right) depicts the properties of stars that are mostly older (younger) than 5 billion years -but see the distribution of the merger times in Fig. 4. Inset histograms show the total stellar mass per galaxy that is formed in-situ in the corresponding periods.
The stars formed before the last major merger (left panel) mostly formed in perfectly circular orbits (circularity ∼ 1). On the other hand, the circularity of their orbits changed as time passed, all the way down to = 0, where they span the whole range of values, with circularity ∼ 0 − 1, i.e. with also non rotationally-supported orbits. This effect is called orbital heating and is manifested in radial and vertical directions. In fact, at = 0, stars can also exhibit negative circularity values, which corresponds to counter-rotating orbits. The heating quantified in the top left panel of Fig. 11 is thought to be due to secular evolution and to more violent orbital changes induced by mergers. In this plot we cannot identify the exact drivers of the orbital alterations nor to assess the timescales when these old stars were heated up to their = 0 levels. Yet, we notice that a tail of the stars in the left panel of Fig. 11, and mostly formed before the last major merger, exhibit circularity values at birth also approaching 0, i.e. random motions: these are at least partially stars that are older than 8 − 10 billion years -and this is qualitatively consistent with the picture suggested by observations and also reproduced by simulations whereby galaxies were dynamically hotter at earlier epochs ( 1.5, Pillepich et al. 2019, and references therein).
Also stars formed during the major merger events and at pericentric passages (top, middle panel of Fig. 11) show broad ranges of circularities at birth, with their circularity distribution at = 0 being even more spread out.
Focusing on the two top right panels of Fig. 11, it is manifest (see also inset histograms) that the stars formed after and during the last major merger are less numerous than the stars formed before it -even just because they could form over shorter periods of time. Yet they are sufficient to secure a disky stellar morphology to the = 0 descendants. Stars formed after the merger (top right) were also born in nearly circular orbits, but, in this case, their circularities barely change during the period between the last major merger and = 0: whereas the latter is comparatively short for disk heating to have a large effect, it should be noticed that for a fraction of the stars in this panel the period elapsed after the last major merger can be as long as ∼ 5 Gyr. Yet, the phenomenology of the top right panel of Fig. 11 quantifies the idea that those MW/M31-like galaxies that experienced a recent major merger are disky at = 0 because of the sustained star formation and because the newly-formed stars are born in circular orbits and have not yet undergone heating.
These results qualitatively agree with the general picture described by Peschken et al. (2020) for = 0 disk galaxies in the Illustris simulation whose last major merger occurred at 1.5: according to that analysis, the stars born before the last major merger form the = 0 spheroidal components and the stars born after the merger constitute a new formed disk.
But how much stellar mass is produced during and after the last recent major mergers? In addition to the insets in the top of Fig. 11, in the bottom left panels we give the fraction of in-situ stellar mass formed in the last 5 Gyr over the total stellar mass (bottom left) and the amount of ex-situ stellar mass brought by the last major merger vs. the in-situ stellar mass in the last 5 Gyr (bottom right). Galaxies with recent major mergers (blue circles) have formed typically a larger fraction of in-situ stellar mass in the last few billion years than their MW/M31 analogues with more ancient last major mergersthis ranges in the median, depending on final mass, between 15 − 30 per cent of the total = 0 mass and the recently formed in-situ stellar masses are, on average, 0.1 − 0.2 dex larger in the recentlymerged population. This indicates that a larger amount of in-situ star formation has indeed occurred because of the recent major merger, Collision angle of the merger orbit vs. mass ratio at max , color coded by the absolute drop in diskyness, D/T, during the merger event, for TNG50 MW/M31-like galaxies that experienced a major merger in the last 5 billion years. The points denote the median of the angles across the 5 to 15 snapshots previous to the merger, whereas the errorbars are the 25th-75th percentiles of such angle distribution). The symbols denote whether the last major merger destroyed (squares) or not (diamonds) the previouslyexisting stellar disk, based on the inspection of the time evolution of the D/T ratios along the main-progenitor branch of each galaxy (see Figs. 8 and 9). The density plots are calculated with a gaussian kernel. Mergers that do not destroy the pre-existing stellar disks tend to populate the parameter space of smaller stellar-mass ratios and larger collisional angles, i.e. spiralling or large impact-parameter orbits (although the stellar mass ratio has a larger impact). and would have not occurred at the same levels in the absence of a recent major merger. At the same time, recent major mergers bring large amounts of ex-situ stars (bottom right panel), a few 10 9 − 10 11 M . This is somewhat necessary, because recent major mergers are more massive in terms of absolute stellar-mass than ancient ones. To summarize, recent major mergers bring both larger amounts of ex-situ stellar mass as well as trigger relatively more in-situ star formation in the last few billion years than in the absence of a major merger. This nicely connects to the average stellar mass growth with time of WM/M31-like galaxies -see Fig. 1, left, and compare black thick vs. blue thick curves (i.e. all TNG50 MW/M31-like galaxies vs. those with a major merger in the last 5 billion years): the latter exhibit a suppressed stellar mass growth at 0.5 than the typical MW/M31-like galaxy, and de facto 'managed' to enter in the stellar mass selection thanks to the mass boost (both in-situ and ex-situ) provided by the recent major merger.

Orbits of the merging galaxies
To understand what determines the disruption or survival of a stellar disk during the merger, we inspect the orbital properties the galaxies follow in their merging process: these can be very diverse. We characterize the orbits of the mergers according to two angles: 1) the angle between the plane of the orbit and the stellar angular momentum of the main progenitor ("orbit plane angle"), and 2) the (acute) angle between the velocity vector of the secondary progenitor with respect to the main galaxy and the position vector between them ("collision angle", as presented in Zeng et al. 2021). These angles are measured at individual snapshots and averaged throughout the 5 last snapshots prior to the merger time (approximately over 800 Myr, to minimize the effect of fluctuations that are sometimes caused by the S halo-finder algorithm). The orbit plane angle allows us to determine whether the secondary galaxy orbit is prograde (angle: 0-90 • ) or retrograde (90-180 • ) with respect to the rotation of the main galaxy. The collision angle determines whether the merging orbit resembles a smooth spiraling (high values) or a head-on collision (low values).
We find (although do not show) that, throughout all the recent major merger events that involve the progenitors of the TNG50 MW/M31-like galaxies, the majority of the mergers that unfold in retrograde orbits (13 from 31 mergers) end up destroying the stellar disks. Also, within our 31 cases, a disk that survives a major merger is more likely to be the result of a prograde than a retrograde orbit. However, these results are tentative and warrant further investigation.
Instead, we can show that it is the combination of orientation of the collision and of stellar mass ratio of the merger that is a good predictor of whether stellar disks are destroyed or not during a recent major merger. This is shown in Fig. 12, where galaxies with recent major merger are color-coded by the drop (in absolute value) of the D/T mass ratio because of the major merger and the latter is characterized by the collision angle and the stellar mass ratio -all of them are major, but the ratios go from 1:4 to 1:1. Galaxies whose stellar disks are destroyed by the major merger and then reform (red squares) are more frequent towards larger stellar mass ratios and lower collision angles, i.e. towards mergers with small impact parameters. On the other hand, the lower the stellar mass ratio, the higher the probability that the disk survives the merger event (see blue diamonds, in general with smaller drops in diskyness) 2 . To separately quantify the effects of the collision angle and of the stellar mass ratio, we perform Anderson-Darling (AD) and Kolmogorov-Smirnov (KS) tests with the null-hypothesis that the two galaxy samples (destroyed and surviving disks) follow the same distribution. For the angles, we cannot reject the null-hypothesis, namely we cannot exclude that the two distributions are in fact indistinguishable (14 (KS) and 18 (AD) per cent significance levels). On the other hand, for the stellar mass ratios, we can reject the null-hypothesis (as their equality is given with very low significance levels, i.e. 1 (KS) and 0.7 (AD) per cent). We conclude that the stellar mass ratios have more significant effects than the collision angle in the survival of stellar disks.

THE PROPERTIES OF MW/M31 ANALOGUES WITH RECENT MAJOR MERGERS
Are the = 0 properties of recently-merged MW/M31-like galaxies different from the other galaxies in the sample, or from the rest of the galaxies in TNG50 in the same mass range? In the following we inspect a selection of global and stellar structural properties of the TNG50 MW/M31-like galaxies and contrast those with recent major mergers to the whole sample. We also juxtapose the measurements from TNG50 galaxies to analogous constraints from observations of the Galaxy and Andromeda: these are indicated in the following plots as magenta and orange areas and bands, respectively, and encompass 2 The reader may notice a galaxy classified as "disk destroyed" but with a very low drop in diskyness. This is Subhalo ID 372754 (top left stamp of  . Mass fraction in kinematically-defined bulges vs. mass fraction in kinematically-defined stellar haloes, for TNG50 MW/M31-like galaxies at = 0. We compare three sub-samples: all MW/M31 analogues (white circles), analogues that experienced at least one major merger throughout their history (light blue), and analogues with recent major mergers (deep blue). Galaxies below the dashed red line have a total stellar mass fraction in both spheroidal components smaller than 0.5 and are therefore diskier. The black dashed line represents an equal fraction of stars in the bulge and the halo. The density plots are calculated with a gaussian kernel. For most TNG50 MW/M31 analogues the stellar halo is more massive than the bulge. Analogues that underwent recent major mergers have average-mass bulges but more massive stellar haloes, compared to the rest of the MW/M31-like galaxies.
the estimates available from the literature (see Table A1 in Pillepich et al in preparation), including systematic differences and measurement errorbars.

Average bulges but more massive and shallower stellar haloes
The origin of stellar bulges in disky galaxies has been debated for at least two decades, with recent observational and theoretical results pointing towards a) a distinction between photometrically-vs. kinematically-defined bulges (e.g. Du et al. 2020, and references therein) and b) an origin of galactic bulges that is not necessarily always connected to mergers, with different pathways for so-called classical and pseudo-bulges (see e.g. Du et al. 2021, and Gargiulo et al. submitted, for TNG50-based results and references therein).
Here we use the kinematic decomposition proposed by Du et al. (2019), based on a Gaussian mixture separation of the stellar particles in the kinematic phase-space, to separate the stellar structural components of TNG50 galaxies in kinematically-defined stellar bulges and stellar haloes. In Fig. 13, we account for all stellar particles that are gravitationally bound and we give the stellar mass fractions in such bulges and stellar haloes, by denoting TNG50 MW/M31 analogues that underwent recent major mergers with blue circles and others with more ancient or no major mergers with light blue and white circles. Galaxies that fall below the red dashed line have a total star fraction in spheroidal components lower than 0.5 and are therefore diskier (106 over 198 MW/M31 analogues, i.e. 54 per cent).
For the entire sample of TNG50 MW/M31-like galaxies the bulge fraction is always (except for one galaxy) below 0.5, and below 0.25 in most cases (∼85 per cent). The mass fractions in stellar haloes cover a wide range, with values also greater than 0.5 (∼ 12 per cent) and reaching values of ∼0.7. When we consider both spheroidal components, there is a non-negligible fraction of TNG50 MW/M31-like galaxies (44/198, ∼22 per cent) where the bulge is more massive than the stellar halo. But, TNG50 returns also bulgeless galaxies (25 galaxies, ∼13 per cent): among these, there are galaxies with major mergers (including recent ones) and without. These latter numbers seem to reasonably agree with previous observational studies of edge-on galaxies (16 per cent found by Kautsch et al. 2006, albeit differently selected). MW-mass galaxies have lower stellar halo mass fractions than M31-mass analogues (median values of 0.27 vs. 0.37), whereas their bulge mass fractions are similar (median of ∼0.15).
For the galaxies with recent majors mergers, the average bulge fraction is comparable to the average values for the rest of the MW/M31 analogues -in disagreement with e.g. the predictions of Puech et al. (2012) of more bulge dominant galaxies with ∼ 0.6 major mergers. A possible interpretative scenario was offered by Hopkins et al. (2010), who pointed to gas richness as a crucial element in limiting bulge growth. However, their stellar halo mass fraction is on average larger, with median halo mass fractions of 0.46 for MW/M31 analogues with recent mergers vs. 0.31 for the complete sample. Again we perform AD and KS to test the null-hypothesis that the two samples (all MW/M31-like galaxies and galaxies having a recent major merger) follow the same distribution, separately for halo and bulge mass fractions. For the halo fraction, the galaxies with recent major mergers have a statistically-distinct stellar halo mass fraction distribution from that of all MW/M31-like galaxies (as the null-hypothesis of equality is given with 0.04 (KS) and 0.1 (AD) per cent significance levels). For the bulge ratios, with 6.9 (KS) and 14.6 (AD) per cent significance levels, we cannot reject the possibility that the two distributions are in fact identical.
The stellar haloes of recently-merged MW/M31-like galaxies are not only more massive, but their stars are also less centrally concentrated, albeit to a somewhat weaker degree. Fig. 14 gives the 3D slopes of the spherically-symmetric stellar profilesfitted with a power law ( ( ) = 0 ) between twice the stellar half-mass radius and 200c of TNG50 MW/M31-like galaxies as a function of the time of their last major merger (top) and of the stellar mass of the secondary galaxy involved in the merger (bottom). Error bars cover one standard deviation errors of the parameter provided by the non-linear least squares fitting routine (python curve_fit). Galaxies with recent major mergers have slightly shallower slopes, with median ([25th-75th] percentiles) of about −4.6 (−[4.3, 5.1]) vs. −4.8 (−[4.4, 5.4]) for the blue and light blue sub samples, respectively. However, this difference is only mildly significant (the significance levels are of 17 per cent and 22 per cent for the KS and AS tests of consistent distributions). This is a distinction that is qualitatively in line but weaker than the findings based on Illustris galaxies by Pillepich et al. (2014). Similarly as found there, TNG50 also predicts the 3D stellar halo slope to depend on galaxy stellar mass: among the MW/M31 analogues, MW-mass galaxies have steeper stellar haloes (−5.0 ± 0.7) than M31-mass ones (−4.5 ± 0.7), a fact that is qualitatively consistent with the estimates of the stellar halo profiles of the Galaxy and Andromeda -see magenta and shaded areas. Furthermore, the more massive the merging companion (within the major merger mass ratios), the shallower the stellar haloes of the descendant galaxy. Error bars are one standard deviation errors of the parameter provided by the non-linear least squares fitting routine (python curve_fit). Magenta and orange shaded area denote current observational constraints in the slopes of the stellar haloes of the MW and M31. In the bottom, we additionally indicate with red squares and black diamonds those recently-merged galaxies whose stellar disk had been destroyed or not by the merger. The stellar haloes of galaxies with recent major mergers have shallower profiles, and among these, the ones whose disk survived the merger but that merged with more massive secondaries are, in general, shallower. We show observational estimates for the MW: for a radius between 25 and 50 kpc, Bell et al. (2008)

Larger fractions and amounts of ex-situ stellar mass
As in the hierarchical growth of structure scenario stellar haloes of MW-like galaxies are mostly formed by accretion (Zolotov et al. 2011;Font et al. 2011;Pillepich et al. 2015), it is not surprising then that MW/M31-like galaxies with recent major mergers also exhibit larger amounts and fractions of ex-situ, i.e. accreted, stars at = 0. This is quantified in Fig. 15, in terms of the distributions of the ex-situ stellar mass fraction across different galaxy samples (top) and the total ex-situ stellar mass as a function of galaxy stellar mass at = 0 (bottom). In both panels, we compare all TNG50 galaxies in At the opposite end, TNG50 MW/M31-like galaxies that never experienced a major merger (white circles) exhibit lower-than average ex-situ stellar masses, given their galaxy mass. the depicted mass range (gray annotations) to TNG50 MW/M31like galaxies (black and light blue lines in the top and white and light blue circles in the bottom) and MW/M31-like galaxies that underwent recent major mergers (blue lines and circles).
The ex-situ fractions range widely: from a few percents to ∼ 65 per cent, with a median value of ∼ 0.19 for TNG50 MW/M31like galaxies. This is in the ball park of the findings by Rodriguez-Gomez et al. (2016) and Pillepich et al. (2018b), with Illustris and TNG100/TNG300 MW-mass galaxies and haloes, respectively. For the sub-sample with a recent major merger, the ex-situ fraction is biased towards the higher end of the distribution, with a median ex-situ mass fraction throughout the galaxy bodies of ∼ 33 per cent. Moreover, at a fixed stellar mass, galaxies with a recent major merger have more ex-situ stellar mass than the average, by 0.3 − 0.5 dex. For the TNG50 MW/M31 analogues with recent major mergers, the total ex-situ mass is always larger than 10 9.8−10 M . Again, the recent major mergers not only trigger star formation (see Figures 10 and 11) but also bring in higher-than average amounts of accreted stars.

Somewhat thicker and hotter stellar disks
Whether = 0 galaxies with thick stellar disks have been born with large stellar velocity dispersions and less flat stellar mass distributions and/or whether their stellar disks puffed up across cosmic epochs, because of secular or externally-triggered processes, is an active field of research (see e.g. Pillepich et al. 2019;Park et al. 2019, for recent works with state-of-the-art cosmological galaxy simulations). Here we do not address this general issue but show that a recent major merger may indeed affect the vertical structure of the stellar disk of a specific selection of galaxies, those with = 0 global properties similar to the Galaxy and Andromeda.
In Fig. 16, we plot the stellar disk heights (top two rows) of all TNG50 MW/M31-like galaxies (light blue circles) and of TNG50 MW/M31-like galaxies that underwent a recent major merger (blue circles). These are shown as a function of stellar disk length (top panels) and of the time of the last major merger (middle panels) for two reasons: a) to control for the influence of galaxy stellar mass and hence disk size; and b) to judge the effects in terms of how long ago a major merger may have perturbed a possibly pre-existing stellar disk. For the disk lengths, for each galaxy projected face on, we fit an exponential profile to the radial distribution of the disk stellar mass surface density in a certain aperture (see Appendix A for details). For the disk heights, for each galaxy projected edge-on, we select only stars in circular orbits in the cylindrical shell between 3.5 and 4.5 times the disk length and fit a double parametric formula to the vertical stellar mass density profile, to allow for both geometrical thin and thick disks. More details about the fitting technique are given in the Appendix A: for each galaxy, we both fit a double hyperbolic secant and a double squared hyperbolic secant profile and represent in the figure the range spanned by both best-fit values. The averages in bins of disk length or time of last major merger (solid thick curves, top and bottom, respectively) run over the mean of the best-fit values from the two functional fitting functions.
According to TNG50, more extended galaxies exhibit thicker stellar disks, particularly for the geometrically thick component. Yet, the galaxy-to-galaxy variation is very large, and in TNG50 we recover example galaxies whose stellar disk structural properties are very similar to both the Galaxy's (magenta shaded areas) and Andromeda's (orange).
From the top panels of Fig. 16, we uncover that galaxies with a recent major merger have, on average, somewhat larger stellar disk heights, at fixed disk length, but this effect is more pronounced for the thick (right) rather than the thin (left) component and for smaller galaxies. From the bottom panels, we also see an increase in disk thickness that mildly correlates with the time of the last major merger, with thicker disks the more recent the last major merger. Although these relationships have significant scatter, they are in place despite the fact that more recent major mergers impact galactic stellar structures that are expected to be thinner and colder (because they exist at lower redshift) than those impacted by merger at higher redshifts, at least at fixed stellar mass (Pillepich et al. 2019). According to AD and KS tests, we find that, for the geometrically-thin component, we cannot reject the possibility that the distributions of the heights of all MW/M31-like galaxies and of those with recent major mergers are indistinguishable (15 (AD) and 20 (KS) per cent significance levels for the null-hypothesis of equality in the case of the double sech 2 ; 3 (AD) and 2 (KS) for the double sech). On the other hand, for the geometrically-thick component, the distributions of disk heights are clearly different (0.7 (AD) and 0.2 (KS) per cent significance levels for the double sech 2 ; 0.3 (AD) and 0.03 (KS) for the double sech).
A similar general picture is in place also when the thickness of the stellar disk is evaluated in a non-parametric way. In Fig. 17, top, we show the stellar half-mass heights of TNG50 MW/M31-like galaxies, measured in the annuli between 3.5 and 4.5 times the disk length: TNG50 MW/M31-like galaxies with recent major mergers have larger stellar half-mass heights than the whole population -and In all panels, vertical colored bars span, for each galaxy, the best-fit values for the scale height obtained from two complementary functional forms: double sech and double sech 2 : light blue vertical bars represent TNG50 MW/M31-like galaxies, blue vertical bars denote those that underwent a recent major merger. Solid thick curves are medians in bins of the quantity on the x-axes of the mean between the best-fit values from the two functional forms. At fixed disk length, galaxies with a recent major merger have somewhat thicker thin and thick stellar disks. Observational estimates for the Galaxy and Andromeda are given as magenta and orange boxes, which encompass, and hence marginalize over, diverse measurements including their errorbars: 1.7 − 2.9 kpc, 175 − 360 pc, and 625 − 1450 pc for the disk length, geometrically thin-and thick-disk heights of the MW (Gould et al. 1996;Ojha 2001;Siegel et al. 2002;Jurić et al. 2008;Rix & Bovy 2013; Bland-Hawthorn & Gerhard 2016) and 4.8 − 6.8 kpc, 900 − 1300 pc, and 2200 − 3400 pc for the disk length, geometrically thin-and thick-disk heights of M31 (Worthey et al. 2005;Barmby et al. 2006;Hammer et al. 2007;Collins et al. 2011). the stellar half-mass height correlates better with the thick-disk height of galaxies than with the thin-disk height. Moreover, when turning to quantify the kinematics of the stellar disks (Fig. 17, bottom), we find that the ratio between the maximum rotational velocity and the vertical velocity dispersion of the stars (as in Pillepich et al. 2019) is lower for MW/M31 analogues with recent major mergers, namely, they are kinematically hotter than the rest of the sample, which concurs at least qualitatively with the structural analysis.
We conclude that MW/M31-like galaxies that have undergone recent major mergers exhibit somewhat thicker and clearly hotter stellar disks, but with a significant scatter. TNG50 does return example galaxies that experienced a recent major merger and have thin and thick disk heights similar to those of Andromeda (orange areas). Interestingly, however, for galaxies with disk length smaller than 3−4 kpc and hence smaller than Andromeda (see magenta vs. orange observational constraints), TNG50 seems to imply that it is improbable, but not impossible, for a galaxy to have emerged from a major merger over the last few billion year and exhibit at = 0 thin-disk heights as low as that of the Galaxy. This is qualitatively in line with the ideas of Toth & Ostriker (1992); Wyse (2001); Hammer et al. (2007), who however excluded the possibility. In fact, such statements strongly depend on how and where the vertical structure of the stellar disks is assessed. For example, in terms of thick-disk or stellar half-mass heights, a few TNG50 galaxies appear as compact and similarly thick as the Galaxy. Moreover, there are two TNG50 galaxies that underwent a recent major merger and whose thin disks has a similar disk length and height as the Galaxy's (Subhalo IDs 528322 and 532301). Finally, when we measure disk heights of TNG50 galaxies at a fixed galactocentric distance of e.g. 8 kpc (irrespective of disk size), we do find a few galaxies that underwent recent major mergers and have a stellar disk as thin as 350 − 400 pc, i.e. similar to the MWs.
These last numbers and considerations indicate that the limited numerical resolution of TNG50 should not affect our scientific conclusions, in that simulated stellar disks can be thinner than the simulation's softening length (288 pc for DM and stellar particles, ≥72 pc for the gas cells -see Section 2). An extended study of the resolution effects on galaxy sizes and heights in TNG50 is presented in Pillepich et al. (2019), where it was shown that, in terms of stellar half-mass heights, stellar disk thickness can be considered to be converged in TNG50 to better than 20-40 per cent. In the top panel, we measure disk heights in a non-parametric way, namely by measuring the stellar halfmass disk heights at 3.5 − 4.5 times the disk length, for comparison to the heights in Fig. 16. Observational constraints on the Galaxy and Andromeda (magenta and yellow areas) are derived from the average relationships in TNG50 between half-mass heights and thin and thick disk heights. In the bottom panel, we show the ratio between the maximum rotational velocity and the vertical velocity dispersion of the stars. Light blue circles represent TNG50 MW/M31-like galaxies, blue circles denote those that underwent a recent major merger, and gray crosses represent the rest of the galaxies in TNG50 in the depicted stellar mass range. Solid curves are medians in bins of the quantity on the x-axes. At fixed stellar mass, galaxies with a recent major merger exhibit thicker stellar disks and lower V/ ratios, i.e. hotter orbits.

Hints of more massive SMBHs, larger gas reservoirs and star formation rates
We conclude this analysis by enunciating (without showing) the differences for a few additional global properties at = 0 of TNG50 MW/M31-like galaxies with and without a recent major merger. Whereas the average stellar-to-halo-mass relation of TNG50 MW/M31-like galaxies with recent major mergers is comparable to the entire sample of MW/M31 analogues, the former seem to be somewhat biased towards larger SMBH masses, at fixed stellar mass. Two effects can explain this: the SMBHs present in the two progenitors were on average larger, even just because they merged at more recent epochs than the others, or the SMBH of the resulting galaxy has accreted more gas and grown more quickly in the post-merger epoch than the SMBHs in the control-sample galaxies. As the effect is not particularly large, we do not investigate further on this here and we refer the interested reader to Aumer, White, Naab & Scannapieco Pil for a discussion on the masses of TNG50 SMBHs vs. that of Sgr A * in the Galaxy.
On the other hand, at fixed galaxy stellar mass, the gas mass fraction and SFR (within 2 1/2, * at = 0) are higher for the galaxies with recent major mergers. This is overall consistent with the picture described above, whereby the recent major merger brought in gas, which in turn cooled down, becoming eligible for star formation. Therefore, not only the gas fractions were already larger at the time of coalescence (Fig. 5), but we also find that MW/M31-like galaxies with recent major mergers are biased high in terms of gas availability with respect to the rest of MW/M31 analogues still at = 0.

SUMMARY AND CONCLUSIONS
We have used the cosmological magneto-hydrodynamical galaxy simulation TNG50 to quantify the mass assembly and merger histories of 198 MW/M31-like galaxies selected at = 0. We have placed a special focus in studying how frequently, and how, MW/M31-like galaxies may undergo a recent major merger (e.g. occurring within the last 5 billion years and with stellar mass ratio > 1:4) and still exhibit a disky stellar morphology at = 0.
The main findings that emerge from our analysis are: • The progenitors of TNG50 MW/M31-like galaxies at ∼ 3 were, on average, ∼40 times less massive in stellar mass than today, but with a large past scatter: 2 dex between the 10th-90th percentiles (7 times less massive in total halo mass, with scatter 0.6 dex; Fig. 1).
• Major mergers are common: 168 of the 198 MW/M31-like galaxies in TNG50 (85 per cent) have undergone at least one major merger throughout their history (or specifically, at 5; Fig. 2).
• TNG50 returns galaxies with stellar mass compatible with the Galaxy and Andromeda and with overall disky stellar morphology at = 0 ( Fig. 3) even in the cases when these have undergone a recent major merger: 31 MW/M31-like galaxies in TNG50 (16 per cent) have experienced at least one major merger in the last 5 Gyr (Fig. 4).
• Galaxies with recent major mergers have interacted with relatively massive companions for significant amounts of times, i.e. on average for ∼1.4 billion years. The companions, i.e. secondary progenitors, are massive objects, with median stellar mass of ∼ 2 × 10 10 M (i.e. average stellar mass ratio of 1:2.5; Fig. 4).
• According to TNG50, there are two main pathways that can lead to a disky MW/M31-mass galaxy at = 0 after a recent major merger: i) the pre-existing stellar disk is destroyed during the interactions and merger with the companion, but reforms (Fig. 8); and ii) somewhat less frequently, the latter does not disrupt the pre-exisisting stellar disk (Fig. 9). Whether the one or the other occurs depends, for example, on the merger configuration, with mergers with larger stellar mass ratios and smaller impact parameters following more frequently the first scenario (Fig. 12).
• In both cases, gas was sufficiently available to either trigger starformation bursts at pericentric passages or at coalescence or both, as well as to sustain prolonged star formation until = 0, with the ensuing (re)formation of a disk of young(er) stars (Fig. 5).
By comparing the = 0 structural and global properties of TNG50 MW/M31 analogues that underwent a recent major merger with those with more ancient last major mergers, we find that the former have, on average: • larger amounts of in-situ stellar mass produced over the last few billion years, namely, as a consequence of the recent merger (Fig. 11); • similarly massive kinematically-defined bulges (Fig. 13); • more massive and somewhat shallower stellar haloes (Fig. 13); • larger amounts and relative fractions of ex-situ, i.e. accreted, stellar mass also at fixed = 0 galaxy mass (Fig. 15); • thicker and hotter stellar disks, but only for the geometrically thick-disk components and for galaxies with smaller disks (Fig. 16); • somewhat more massive SMBHs at fixed = 0 galaxy mass, and larger gas mass reservoirs and higher star formation rates even at = 0.
These results suggest that it may be possible to associate a probability for an observed galaxy to have experienced a recent major merger on the basis of some key structural and global observable properties (Zhu et al. 2021;Eisert et al. 2022).
Importantly, the results quantified in this paper align well with the current observational constraints on the properties and on the recent past assembly histories of the Galaxy and Andromeda. As anticipated in the Introduction, whereas the last major merger of the Milky Way may have occurred as early as 8 − 11 billion years ago, Andromeda's last major merger may have happened only a few billion years ago. As indicated throughout the paper by reporting known observational properties from the literature, the Milky Way exhibits a steeper stellar halo profile, a thinner and colder stellar disk, and a higher star formation rate -given its mass -than Andromeda: all this is consistent with the average phenomenology of TNG50 MW/M31-like galaxies with ancient and recent last major mergers.
With this paper we clearly signal that, with current state-of-theart cosmological galaxy simulations that encompass wide ranges of merger histories within the ΛCDM hierarchical growth of structure, there is no need to simulate galaxies with quiet recent merger histories to obtain galactic stellar disks at = 0. This approach had been frequently adopted over the last decade with zoom-in cosmological simulations of MW-or M31-mass galaxies, whose parent DM haloes had to be chosen from lower-resolution, DM-only volumes for resimulation. Whereas this choice has almost always been dictated by the difficulties of the overcooling problem and by practical strategies (quieter merger histories imply faster computing times), it has the drawback that, in practice, no past and recent zoom-in simulations could really address the question of whether the absence of a major merger since 1 is a necessary condition for a galaxy to have a stellar disk as thin as the Galaxy's (i.e. 175 − 360 pc and 625 − 1450 pc for geometrically thin and thick components, respectively). Additionally, high-resolution galaxy simulations that are suitable for understanding the formation and evolution of Andromeda have been, by imposition or necessity, exceedingly rare, at least until recently and until TNG50. In fact, TNG50 produces a few example galaxies whose more detailed stellar disk structures are also compatible with the Galaxy's and Andromeda's. Even more interesting, TNG50 produces example galaxies that did experience a recent major merger and have thin and thick disk heights similar to those of Andromeda (0.9 − 1.3 and 2.2 − 3.4 kpc, respectively). Moreover, according to TNG50, it does seem improbable, but not impossible, for a galaxy as small (in disk size) as the Galaxy to have emerged from a major merger over the last 5 billion years with thin-disk height as small as the Galaxy's. However, these statements do depend on where within the disk and how the vertical structure is assessed.
We conclude this discussion by highlighting that important structural changes and star formation episodes can be triggered in the progenitors of MW/M31-like galaxies also prior to coalescence with a major companion, i.e. during pericentric passages (see Figs. 6 and 10). It is hence to be expected, and to be searched for, that the Milky Way's disk and stellar halo may exhibit perturbations (such as waves and vertical spirals) that are linked to the passage of the Large Magellanic Cloud, as recently pointed out with observational data by e.g. Vasiliev et al. (2021) and Conroy et al. (2021). Our findings also offer a qualitative glimpse into what our Galaxy may experience in just a a couple of billion years in the future (Cautun et al. 2019), when the Large Magellanic Cloud will eventually merge with it.