DOES M31 RESULT FROM AN ANCIENT MAJOR MERGER?

Our nearest giant neighbor, M31, has attracted a considerable amount of interest since the discovery of many large-scale structures that surround its outskirts. These prominent structures include the Giant Stream (infall of 1.5 × 10⁸ M_☉ stellar mass; Ibata et al. 2001) and the gigantic thick disk containing about 10% of the disk luminosity (Ibata et al. 2005). The halo of M31 is haunted by not less than 16 substructures that led Tanaka et al. (2010) to conjecture that they are due to as many accretion events involving dwarf satellites with mass 10⁷–10⁹ M_☉ since z ∼ 1. Observations at radio wavelengths (Westmeier et al. 2005; Thilker et al. 2004, and references therein) reveal high velocity clouds, whose location and kinematics partly follow the Giant Stream. Moreover, a 10 kpc pseudoring inserted in the M31 disk dominates star formation (Baade & Arp 1964), H i gas (Roberts 1966), molecular gas (Nieten et al. 2006), and dust emission (Gordon et al. 2006).

Simulations of the numerous structures in the M31 outskirts have always assumed them to be caused by minor satellites. M32 is understood to be the perturber of the spiral arms (Byrd 1978, 1983), while NGC 205 has been modeled as the cause of the warp seen in the optical and H i disks of M31 (Sato & Sawa 1986). The thick disk has been modeled by earlier disruption, several Gyr ago, of dwarf galaxies on prograde orbits that are coplanar with the disk (Peñarrubia et al. 2006). Simulations of the Giant Stream assumed a collision with an unknown satellite ∼0.7 Gyr ago (Font et al. 2006; Fardal et al. 2008). The 10 kpc ring could also be reproduced by an interaction with M32, assuming a polar orbit (Block et al. 2006).

In spite of their success in reproducing M31 structures, several of these simulations may be speculative because the M31 satellite orbits are currently unknown (e.g., Fardal et al. 2009). More problematic is the fact that stars in the Giant Stream (Brown et al. 2007) have ages from 5.5 to 13 Gyr, which is difficult to reconcile with a recent collision (e.g., Font et al. 2008). If the substructures are formed from different progenitors, why do they show obvious similarities in metallicity (e.g., Ferguson et al. 2005)? The main motivation for only investigating minor mergers is to preserve the M31 disk age (Mori & Rich 2008), but the M31 disk is not necessarily very old and permanent, as its associated stellar clusters display young to moderate ages (< 5–7 Gyr; Beasley et al. 2004), and most of the stars in the outskirt structures are older than this. One should not discard the possibility that the gigantic structures may have been formed at the same time or even earlier than the disk. In fact, many co-workers in the field (Rich 2004; Ibata et al. 2004; Brown et al. 2006) have hypothesized a possible major merger in the past history of M31.³ More recently, the kinematics of the M31 globular system has been attributed to an ancient major merger (e.g., Bekki 2010).

Understanding the nature of M31 certainly has an impact in cosmology because it and the Milky Way are the only two giant spirals in our immediate neighborhood. In contrast to the Milky Way, M31 has properties (absolute K luminosity, angular momentum, [Fe/H] in the outskirts) similar to the average of large spirals having the same rotation velocity (Hammer et al. 2007). Half of the progenitors of spiral galaxies in this range were not relaxed 6 billion years ago, according to the detailed studies of their morphologies (Hammer et al. 2005; Neichel et al. 2008; Delgado-Serrano et al. 2010) and kinematics (Puech et al. 2008; Yang et al. 2008). Hammer et al. (2009) verified that various merger phases can reproduce quite well the observed morphologies and kinematics of these non-relaxed galaxies. This agrees with the disk rebuilding scenario (Hammer et al. 2005), according to which most spiral disks have been rebuilt following a major, gas-rich merger during the past 8–9 Gyr. High gas fractions have been shown to be essential during this process (Hopkins et al. 2008, 2009b). Understanding the different substructures of M31 as resulting from a single event such as a gas-rich major merger is an important step in validating or dismissing such a channel for spiral disk formation.

The goal of this study is to investigate whether a past major merger, instead of multiple minor mergers, can reproduce most of the peculiarities of M31. The Λ cold dark matter cosmology (H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7) is adopted throughout the paper.

In the following we assume that M31 results from a past major merger. Our goal is to simultaneously reproduce (1) the M31 disk at its rotation velocity, and the bulge, with B/T = 0.28, (2) the presence of a very extended, co-rotating thick disk, (3) the 10 kpc dust, H i, star-forming ring, (4) the Giant Stream, and (5) the observed age distributions of stars in the different substructures. Such a choice may appear somewhat subjective, although it has been adopted after thoroughly reviewing many of the M31 properties. First, because the M31 system is very complex to model, we need to limit the number of features that can be reproduced. Second, the limitations of our model (number of particles) prevent us from describing very small structures, for example, the double nucleus of M31 (but see Hopkins & Quataert 2010). Third, we aim at reproducing the zero- and first-order features of each substructure of M31, including their stellar age, mass, and kinematics (if they are available from observations), but not the detailed morphology of each of them (see discussion, Section 5).

According to Hopkins et al. (2009b, see their Figure 7), a mass ratio from 0.3 to 0.5 is required to reform a new bulge with B/T ∼ 0.3 in the remnant galaxy for gas fractions⁴ ranging between 0% and 50%. However, reforming a significant thin disk requires that enough gas be preserved after the fusion of the cores: it is mostly made of stars from the gaseous thin disk that immediately reforms after the collision (Barnes 2002; Abadi et al. 2003; Governato et al. 2007; Hopkins et al. 2009a). Thus, reforming an Sb galaxy like M31 probably requires a gas fraction in excess of 50% in the progenitors. The total baryonic mass of M31 is 1.1 × 10¹¹ M_☉ (Hammer et al. 2007), and at first approximation we can assume that this is the value for the sum of the two progenitor masses. In order to prevent a too violent relaxation at the center of the main progenitor, we adopt a prograde–retrograde orientation of the spin axes, which has been shown to be more favorable in rebuilding a thin disk after a major merger (Hopkins et al. 2008).

The presence of a gaseous ring favors a polar orbit. As shown by Mori & Rich (2008), the Giant Stream may result from the returning material of the tidal tail formed just before the last passage of the satellite. A similar although longer lived tidal tail is predicted for a 3:1 major merger with a polar orbit, which is associated with the passage of the secondary galaxy just before fusion. During the first passage and until the fusion, star formation is especially enhanced in the secondary galaxy that is harassed by the main one (e.g., Cox et al. 2008). If the Giant Stream was associated with returning particles from the tidal tail formed during the second passage (near fusion), it would contain mostly stars older than the epoch of the fusion; this is because the star formation cannot hold for a long time within tidal tails due to their expected dilution.

We thus propose from Figure 1 a chronological history of the different structures in M31, as this figure can be used as a clock for determining the occurrence of merger phases. The star formation history of the whole system is enhanced during the first passage until the fusion and then at the fusion itself (see Cox et al. 2008). During a merger event, most of the gas and stars in the remnant outskirts have been deposited by tidal tails formed during the first passage and later during the fusion of the cores. A few hundred billion years after its formation, the tidal tail dilutes, provoking a natural quenching of the residual star formation (see Wetzstein et al. 2007). Thus, the age of the material brought by tidal tails provides, with a relatively small delay, the date of both the first passage and fusion times. In the following, we assume that the first passage occurred from 8.5 to 9 Gyr ago, and that the corresponding tidal tails are responsible for the halo enrichment seen in the 21 and 35 kpc fields, without significant star formation more recent than 8.5 Gyr. The thick disk has a star formation history comparable to that of the Giant Stream and is also generated by material returning to the galaxy mostly from tidal tails generated at the fusion. Because their youngest significant population of stars has ages of 5.5 Gyr, the delay between the first passage and fusion ranges between 3 and 3.5 Gyr. This could be accommodated for by relatively large impact parameters (20–30 kpc).

Zoom In Zoom Out Reset image size

From the above, we can settle the initial conditions for a major merger assumed to be at the origin of M31. Table 1 describes the adopted physical parameters of such an interaction. Given the huge amount of observational data, it is beyond the scope of this study to reproduce the details of all the numerous structures in the M31 system. Instead, our aim in using hydrodynamic simulations is to determine whether or not these numerous substructures can be attributed to a single major merger in the past history of M31.

We used the GADGET2 hydrodynamical code (Springel 2005) and initial conditions similar to that of Cox et al. (2006, 2008). For the dark matter, we adopted a core density profile as in Barnes (2002), with a core size of 5.3 and 3.06 kpc for the main galaxy and the satellite, respectively. We verified that our results are not significantly affected by changing the density profile to a Hernquist model that is quite similar in central regions to a Navarro–Frenk–White model (see Section 5). Concerning the implementation of feedback and cooling, we have verified step by step our ability to reproduce both the isolated and merger cases for all the different combinations of feedback and cooling in Cox et al. (2008), and as such, all parameters used in our simulations are very similar to those of Cox et al. (see J. L. Wang 2010, in preparation). Usually, the free parameters describing the star formation efficiency and the feedback strength are chosen in order to match the Schmidt–Kennicutt law between star formation and gas surface densities (see Cox et al. 2006). Given the relatively large scatter of this relation, a large number of combinations can be accommodated (Cox et al. 2006, 2008). In the following, we explored several combinations that are in agreement with this relation, but with the additional constraint of preserving a significant gas reservoir in the progenitors, before fusion (see Section 4.2). Properties of the progenitors are listed in Table 1, and these galaxies have been generated to follow the baryonic Tully Fisher relation (see Puech et al. 2010, for an argument in favor of a non-evolving relation). Their atomic gas content has been assumed to be 3 times more extended than the stars, as also adopted by Cox et al. (2006); see also an observational support from van der Kruit (2007) and references therein. We have to preserve enough gas immediately after the fusion to preserve the formation of a significant thin disk with more than 50% of the baryonic mass.

During our investigations to optimize the orbital parameters, we run the GADGET2 code with 150,000 particles. In most simulations, we have ensured that the mass of dark matter particles is no more than twice that of initial baryon particles, i.e., gas or stars in the progenitors. During star formation events controlled by the gas volume density (and following the Kennicutt–Schmidt law), we have also limited the number of newly formed stars per gas particle to 3. This is an important limitation to the simulation as we have verified that when newly born stars are relatively too small in mass, they are artificially scattered due to their encountering heavy dark matter particles. Similarly, we need to keep the mass of dark matter particles low enough to avoid non-physical disruption of the newly reformed disk at the end of the simulation. In the following, we have tested that the decomposition of the newly formed galaxy does not depend on the adopted number of particles within the range of 150,000–800,000.

However, simulating the first-order properties of faint structures like the Giant Stream requires a high particle number: the Giant Stream is made of 3 × 10⁸ M_☉, which corresponds to 1/366th of the total baryonic mass of the galaxy, and allows only a few hundreds of particles in the Giant Stream for a 150,000 particles simulation. Having fixed the range of parameters to reproduce the observed decomposition of M31 in the bulge and the thin and thick disks, we used a limited number of simulations with 300,000–800,000 particles to better identify the formation of faint structures.

Table 1 lists the adopted range of parameters which have been investigated (Column 2) and adopted (Column 4) for several models considered in this study. These models are all part of the same family of orbital parameters, and they differ essentially by various adjustments of the star formation history or small variations of orbital parameters. We have begun with a model with moderate gas fractions (40% and 60% in Gal1 and Gal2, respectively) to verify whether a significant thin disk can be formed after fusion, assuming somewhat extreme conditions (e.g., high feedback and very low star formation efficiency) to prevent star formation before the fusion epoch. We then realized that most of the gas preserved before fusion forms a very thin disk (see also Abadi et al. 2003), which also hosts some of the pre-existing stars before fusion. However, the thin disk fraction is generally not much larger than the available gas fraction before fusion. Because about 65% of M31 stars lie in the thin disk, we have generated for our current models a higher gas fraction that allows us to test less extreme conditions for feedback and star formation efficiency. Indeed, such gas fractions are quite common at z = 1.4 (see Daddi et al. 2010), a redshift which corresponds to 9 Gyr ago, just before the interaction between the M31 progenitors (see Section 2). Finally, the models adopted in Table 1 (Column 4) are developed to optimize the reproduction of the decomposition of the galaxy in the bulge and the thin and thick disks, and the disk scale length (see Sections 4.1 and 4.2), the presence of the 10 kpc ring (see Section 4.3), the Giant Stream and its kinematics (see Section 4.4), as well as the fraction of stellar ages in most of these substructures (see Sections 4.2 and 4.4). An overview of the goodness of each model is provided in Section 5.

Figure 2 shows the different steps of the merger for one of our models. At the first passage (at T = 0.7 Gyr), tidal tails are formed and enrich the halo, although they are mostly diluted at the fusion epoch. Part of the material ejected during the second passage and the fusion progressively return to the galaxy after the fusion. The material associated with the tidal tail pointing to the bottom of the galaxy (see panel at T = 4.2 Gyr) possesses an angular momentum very different from that of the newly formed disk. This perturbation may explain such faint features like the Giant Stream, and could be persistent for several Gyr after fusion (see Section 4.4).

Zoom In Zoom Out Reset image size

4.1. Decomposition of the Newly Formed Galaxy in Sub-components

We choose to decompose the newly formed galaxy in three components, according to the relative strength of the angular momentum in the direction of that of the thin disk, which is approximately the orbital angular momentum. A very accurate determination of the thin disk angular momentum is possible because most of the preserved gas particles before the fusion naturally form a thin disk (see also Abadi et al. 2003) in which most of the star formation occurs after fusion. We have used the stars born 0.5 Gyr after the fusion to determine the angular momentum of the newly formed disk (see Figure 3). With this technique, we obtain an accuracy of ±1° for determining both the disk position angle (P.A.) and inclination.

Zoom In Zoom Out Reset image size

Figure 4 (top panels) shows the ratio of the angular momentum along the y-axis to the total angular momentum (J_Y/J_total) as a function of the x-axis along which the thin disk is projected (see Figure 3). The final disk rotates anti-clockwise (see Chemin et al. 2009), and the thin disk (illustrated in Figure 3) is easily recognizable in Figure 4 top left panel for values ranging from −0.9 and −1, which correspond to the gaseous thin disk. Other recently formed stars are located within a compact structure without preferred angular momentum orientations, which corresponds to the bulge. Ibata et al. (2007) pointed out the confusion in the literature between the spheroid component and the thick disk. They revealed that the minor-axis region between projected radii of 7 kpc < R < 18 kpc is strongly affected by a rotationally supported thick disk that corresponds to approximately 1/10th of the thin disk stellar mass.⁵ For a better comparison with observational data, we assume that the thick disk is made by all particles (except those of the thin disk) whose angular momentum is dominated by rotational motions along the thin disk. It means that J_Y has to be larger than the combination of other components of the angular momentum (i.e., $\sqrt{J_X^2+J_Z^2}$), i.e., J_Y/J_total should be larger than −1/$\sqrt{2}$ = −0.707, as illustrated by the dot-dashed lines in Figure 4. We have generated the decomposition proposed by Abadi et al. (2003) for comparison. In their scheme, the thick disk is the residue after subtracting the bulge and the thin disk, assuming a bulge symmetrically distributed around J_Y = 0. While the thin disk fraction is found to be very similar with both approaches, the residual thick disk from Abadi et al.'s method is almost twice that provided by our method. We notice that a significant fraction (up to 50%) of such a residual thick disk is made of particles which are not dominated by rotation around the y-axis, which well explains the discrepancy. For comparison with the Ibata et al. (2005) observations of a rotationally supported thick disk, we keep our decomposition as described above.

Zoom In Zoom Out Reset image size

4.2. Thin and Thick Disks: Feedback Prescriptions

Table 2 describes how the decomposition depends on the choice of orbital parameters. From the second to the fourth sets of parameters, only the mass ratio varies and the bulge fraction increases as the mass ratio decreases as shown by Hopkins et al. (2009a). However, we have also tested a higher mass ratio (first set of parameters) with a different set of inclinations for the progenitors. Thus, the bulge fraction may also depend on other orbital parameters, here the relative inclination of the progenitors. For example, it reaches a maximum when the difference of inclination between the two progenitors is close to 180°, possibly because more stars take radial orbits and fall into the bulge. It will be shown later in Section 4.4 that this resonance has a more considerable effect on the matter that is stripped from the satellite to the tidal tail at the second passage.

Table 2. Decomposition in Mass (Unit = 10¹⁰M_☉) of the Newly Formed Galaxy 9 Gyr After the Simulation, for Various Parameters of Our Modeling

Parameters	Comp.	Thin	Bulge	Thick
		Disk		Disk
rpericenter = 24	Stars	2.66	2.51	0.85
Mass ratio = 3.5	Gas	3.24	0.15	0.09
Gal inc = 65 and −110	Fraction	62%	28%	10%
rpericenter = 24.8	Stars	2.72	2.31	0.78
Mass ratio = 3.0	Gas	3.34	0.12	0.11
Gal inc = 65 and −89	Fraction	65%	26%	9%
rpericenter = 24.8	Stars	2.38	2.67	0.95
Mass ratio = 2.8	Gas	2.95	0.13	0.10
Gal inc = 65 and −89	Fraction	60%	29%	9%
rpericenter = 24	Stars	2.25	2.84	0.88
Mass ratio = 2.5	Gas	2.91	0.12	0.09
Gal inc = 65 and −89	Fraction	57%	32%	10%
Same as above	Stars	3.62	2.61	1.22
Feedback	Gas	1.32	0.07	0.12
Changed at 3.5 Gyr	Fraction	55%	29%	14%

Notes. For all models but the last one, high feedback has been maintained during all the simulations. For the last simulation, feedback has been abruptly dropped from high to low values (from Cox et al. 2008), 3.5 Gyr after the beginning of the simulation and ∼0.7 Gyr after the fusion.

Download table as:  ASCIITypeset image

Almost all the simulations we have generated are re-forming significant thin disks containing more than half of the baryonic mass. This is due to our feedback prescriptions that always preserve a significant fraction of gas just before the fusion. As stated by Cox et al. (2008), very little is known about the requisite conditions enabling the star formation to occur. An important prerequisite, however, is the presence of stable gas-rich galactic disks at high redshifts, that could be the progenitors of M31. Cox et al. (2008) have tested feedback for an isolated disk galaxy with an initial gas fraction of f_gas = 0.52. They found that the gas fraction is mostly preserved within 1 Gyr if median to high feedback prescriptions are adopted, while low feedback decreases the gas fraction to 0.37 in the same amount of time.

In our current models, we have assumed f_gas ∼ 0.65: the re-formation of a significant thin disk implies that most of the gas is not transformed into stars before fusion. This calls for relatively high feedback, at least intermediate between the median and high values of Cox et al. (2008). We have verified that for half the value of high feedback of Cox et al. (2008), most of the initial gas is preserved before the fusion. However, re-forming a significant stellar thin disk requires a significant change of the feedback at or after the fusion, in order to transform stars from the gas in the newly reformed disk. It results that from the merger of realistically gas-rich galaxies at z ∼ 1.4, we may reform a galaxy resembling M31, with the condition of a transitory history of the feedback during the merger, from high to low values. There could be some theoretical grounds favoring such a history, however. For example, in a gaseous-rich and almost pristine medium, first supernovae are likely massive and generate high feedback as they may delay star formation for up to 100 Myr (e.g., Bromm et al. 2009). The fusion of two gaseous-rich galaxies corresponds to a severe mixing of most of their components, accompanied by large star formation rates. Then, the metal abundance in the remnant becomes larger and more homogeneously distributed than it was in the progenitors before fusion. The occurrence of very massive, primordial supernova is unlikely in such a mixed medium leading to a possible transition from high to low feedback. Before the fusion there are many fewer exchanges between the two interlopers, and it might be realistic to assume a negligible change in the feedback history.

Of course, the above is mostly made of conjectures in the absence of observations of gas-rich galaxies at high redshift showing no or very small amounts of star formation, i.e., the high-redshift counterparts of the present-day low-surface brightness galaxies. We may, however, verify whether our adopted star formation history is consistent or not with the star formation history revealed in each M31 substructure displayed in Figure 1. The thick disk is indeed mostly made by matter returning from the tidal tails. It is comprised of approximately 80% stars older than 8 Gyr, the ages of the other stars being mostly from 5 to 8 Gyr. Such a distribution is well reproduced by our model with high feedback values (86% of old stars) or with half this value (80%). Because the thick disk shows a similar fraction of stars in the two bins set by Brown et al. (2008) at 5.5 and 7 Gyr, it is probable that there are no significant changes of the star formation history between the two corresponding epochs, i.e., between the first and second passages.

Within these prescriptions, our models predict quite well the fraction of each sub-component of M31 (see Table 2). Note also that the modeled thin disks show rotation curves and scale lengths that are in good agreement with observations (see Figure 5). In the following sections, we examine whether this model could also reproduce other large-scale structures such as the ring and the Giant Stream.

Zoom In Zoom Out Reset image size

4.3. The 10 kpc Ring

The decomposition of the light profile was done by fitting the radial profile with a three-component model (bulge and thin and thick disks). This provides a thin disk with a Sérsic index of 1 and a scale length varying from 4.8 to 5.1 kpc for models with mass ratio ranging from 3 to 2.8. Such values are slightly smaller than the 5.8 ± 0.4 kpc that was adopted by Hammer et al. (2007) for the M31 disk. All these decompositions (see Figure 5) reveal the presence of a prominent ring at 10 kpc which corresponds well to the observed ring. Similarly, all models reproduce the 10 kpc gaseous ring (see Figure 6) as seen in H i observations. This is not unexpected as all the orbits are close to polar which favors the formation of such prominent and persistent structures. In fact, H i observations do not provide many constraints on our modeling. They only discard parameters for which the disk is too warped, i.e., implying that the difference between the two inclination angles of the progenitors should not be too close to 180°. Moreover, we notice (see Figure 6) that with a constant high feedback during the simulation, there is still gas within the bulge of the remnant galaxy. However, this gas is consumed in models that assume a transition to low feedback immediately after the fusion.

Zoom In Zoom Out Reset image size

4.4. The Giant Stream and Its Kinematics

The Giant Stream is a very faint stellar structure, which extends from 50 to 80–100 kpc from the M31 center. We have modeled it as being caused by particles coming back from the tidal tail formed just before fusion. Figure 2 describes the formation of such a structure that is aligned along the trajectory of the satellite which falls into the mass center at the fusion, 4.5 Gyr after the beginning of the simulation. We verified that, within the family of orbits we choose, the strength of the tidal tail (see Figure 2, panel at 4.2 Gyr) depends on the inclination of the progenitors relative to the orbital angular momentum (see Table 1). Optimal values are found for large values of the inclination of the main galaxy (>50°) and especially for differences between the inclination of the two progenitors between 140° and 170°, i.e., not too far from the resonance at 180°.

Modeling the Giant Stream is not an easy task as the structure is very faint compared to other large substructures of M31, even the thick disk. As in Section 4.1, we make use of angular momentum properties to identify peculiar structures able to persist after the remnant phase (see Figure 7).

Zoom In Zoom Out Reset image size

Figure 7 shows the distribution of the angular momentum for the thick disk (top left panel) and at the galaxy outskirts, 5 Gyr after fusion. The top right panel shows two peaks which correspond to the two tidal tails formed after the second passage. The tidal tail associated with the main galaxy (see Figure 2; at 4.2 or at 4.5 Gyr) is spiraling around the thick disk and its angular momentum peaks at Φ = 290°, which is indeed part of the thick disk (see top left panel). This peak disappears in the bottom left panel for which the thin and thick disks have been removed. The residual particles form another well-identified peak at Φ = 225°. We have verified that they correspond to particles returning from the tidal tail of the satellite: in the bottom right panel, selected particles of the tidal tail of the satellite indeed show a very similar distribution with a peak at Φ = 225°, which is offset by 45° from the disk angular momentum. The selection of these particles has been done in a rather crude way (see Figure 2; panel at T = 4.2 Gyr): we simply "cut the tidal tail" at the time of its formation, up to the edge of the satellite. However, this crude selection misses some particles, for example, those which are stripped from the satellite at later phases (e.g., third passage, see Figure 2; panel at 4.5 Gyr). Comparing the counts of particles having an angular momentum around the Φ = 225° peak, we find that our selection recovers 90%, 70%, 40%, and 15% within shells of 100–170, 50–100, 30–50, and 20–30 kpc, respectively. As expected, particles at the very edge of the tidal tail are the most difficult to pre-select. Five gigayears after fusion those particles have returned to trajectories close to the galaxy center.

Tracing the trajectories of selected particles in the tidal tail is very instructive: they correspond to an important component of the angular momentum which is at odds with the rotation of both the thick and thin disks (see Figure 8).

Zoom In Zoom Out Reset image size

Figure 8 describes the temporal evolution of the trajectories of tidal tail particles returning to the galaxy. When falling to the galaxy center, the particles are looping around the galaxy with a pericenter at ∼25 kpc. These particles are describing several elliptical orbits—or loops—around the galaxy potential. We identified four of these loops, although they might be more numerous. Figure 8 shows that with increasing time, they show a precession around the galaxy as they grow. This is because stars coming back from the tidal tail are returning from higher elevation at a later epoch and are coming back with increasing energy. We also note that the loops are somewhat thick and, as such, they could better be called tores.

The permanent rain of stars from the tidal tail ensures very long-term streams of stars within the galaxy outskirts well after the remnant phase. Infalling stars are trapped into elliptical orbits with the remnant galaxy right at the foci, which make the orbit stable. These orbits are those expected for extremely small satellites (single stars!), in the absence of tidally induced forces. Loops are expected to persist for several billion years⁶ after fusion and, moreover, they are permanently fed by new particles coming from the tidal tail. These streams and loops are all within a common plane, 45° from the P.A. of the rejuvenated disk as is evidenced by Figure 8 (panels (h) and (i)). This is our proposed model for the formation of the Giant Stream, which indeed points toward the center of the galaxy. Let us now examine how this model is compliant with several observations of the Giant Stream properties.

The kinematics of stars follows the trend evidenced by Ibata et al. (2004) and Gilbert et al. (2009): stellar particles in the streams reach large negative values of their heliocentric velocities toward the center (see panel (i) of Figure 8). Figure 9 shows the distribution of heliocentric velocities for stellar particles selected in the region of the Giant Stream within Φ_obs from −90° to −140°. Observations (see Figure 1 of Ibata et al. 2004 and Figure 5 of Gilbert et al. 2009) show a broadening of the velocity distribution toward the galaxy center that is well reproduced by our simulations. Besides this, our model also well reproduces the distribution of stellar ages observed by Brown et al. (2008), as we find 80%, 20%, and 0% of stars older than 8 Gyr, from 5 to 8 Gyr, and younger than 5 Gyr, respectively. These fractions can be regulated by changing the feedback prescriptions (see Section 4.2).

Zoom In Zoom Out Reset image size

We have attempted to reproduce the photometric properties of the Giant Stream. The difficulty of such an attempt is due to the fact that the Giant Stream is very faint when compared to other sub-structures of M31, which is a problem also encountered by observers. Even with 500,000 particles, the number of stellar particles in the Giant Stream is below 1000 which obviously limits a detailed reproduction of the structure. Within an angle of Φ_obs from −90° to −140° and from 25 to 80 kpc, the total mass of the Giant Stream is 4.8 and 4.4 10⁸ M_☉ in the two examples illustrated in Figure 10. Using other variations of the model parameters, we find a range from 2 to 5 × 10⁸ M_☉, which easily includes the observational value of 3 × 10⁸ M_☉ (see Mori & Rich 2008).

Zoom In Zoom Out Reset image size

All our models (Column 4 of Table 1) show streams of stars which belong to a plane that intersects the observational plane in a direction 45° from the thin disk P.A. We verify that to optimize the reproduction of the observed features in position and strength, we can use two different kinds of symmetries. First, we can rotate the whole system around an axis perpendicular to the thin disk, as this would not affect the galaxy decomposition in sub-substructures, nor the kinematics of thin and thick disks. Second, we find that gravity is not sensitive to a complete inversion of the initial (x, y, z) coordinate system. This has the advantage of being able to reproduce the Giant Stream with many variations of our model parameters, but it also has the inconvenience of providing a very large amount of space parameters to investigate. Figure 10 illustrates some examples of this exercise for two models with 500 K particles for which the Giant Stream mass is slightly larger than the observed one. Much larger numbers of particles are certainly required to simultaneously reproduce all the morphological details of the thick disk, although we generally reproduce the NE and W shelves quite well.

By construction, all models reproduce the distribution of stellar ages in the various substructures quite well, as they are described in Figure 1. This is especially true for both the thick disk and the Giant Stream (see Sections 4.2 and 4.4) for which our predicted compositions of stellar ages are from 10% to 20% of intermediate age stars (5–8 Gyr), the rest essentially being made of older stars. However, we experienced some difficulties in reproducing the distribution of the 21 and 35 kpc "halo" fields of Brown et al. (2008). We find a fraction of intermediate-age stars ranging from 5% to 15%, while according to Figure 1, it should be near or below 5%. It would be tempting to use this property for selecting the best model of M31. However, in our models intermediate-age stars are related to low-mass particles, three of them being produced by a single gas particle. It may lead to an artificial scattering of those intermediate-age star particles (see Section 3), which could be only solved by simulations with a much larger number of particles than 500 K.

All models as defined in Table 1 (Column 4) can reproduce with quite good accuracy the bulge and thin and thick disk fractions, as well as the ring and the shape, the mass, and the kinematics of the Giant Stream. The examples shown in Figure 10 are cases for which the Giant Stream stellar content is slightly larger than the observed one, for illustrative purposes. We also note that the shapes of the predicted thick disk and Giant Stream can significantly change—within a given set of orbital parameters (see Table 1)—either because they can be modified through a rotation around the disk axis or because they could fluctuate with the assumed time after the beginning of the simulation. It is beyond the scope of this paper to search for the best model that reproduces all the details of M31 substructures. This would require launching several hundreds of simulations with several millions of particles. We have, however, learned the following from this analysis.

• 1.  
  A thin disk such as that of M31 can be reproduced with a combination of ∼65% of gas in the progenitors and a star formation history that prevents gas consumption before the fusion; low efficiency of star formation is thus predicted in primordial gaseous phases, while when the medium has been sufficiently enriched and mixed during fusion, the gas is more easily consumed; this transition has essentially been modulated through the assumed feedback history.
• 2.  
  The same star formation history also reproduces a rotationally supported thick disk with a stellar composition similar to what is observed.
• 3.  
  Fractions of bulge/thin disks and thick disks similar to what is observed within 10% accuracy.
• 4.  
  These polar orbits always reproduce the observed 10 kpc ring, and in most cases, the gas map is very similar to the observed one (see Figure 6), especially when the gas in the galaxy center is allowed to form stars after the galaxy fusion.
• 5.  
  The predicted Giant Stream is made of stellar streams due to returning stars from a tidal tail formed prior to fusion from stripped particles of the satellite; these stars are drawing streams along loops which are inserted into a thick plane that is inclined 45° relative to the disk P.A. (see Figure 8). Such streams are very long-lived because returning particles are trapped into such loops for periods larger than several gigayears, and also because the streams are permanently fed by new particles falling back from the tidal tail at all epochs.

The success of modeling M31 properties as resulting from a past major merger event is patent. Its predictive power is, however, limited by the complexity of the system as well as the very demanding number of particles required to reproduce in detail all substructures of M31. An important limitation could also be the accuracy of the calculations, which requires tracking positions and velocities of particles with an unprecedented accuracy for almost a Hubble time. In fact, we cannot claim to have determined the best major merger model of M31, although it is likely within the limits discussed below.

We also did not try to investigate the whole parameter space to which the above results hold, although examination of Table 1 led us to suspect that it is likely large. However, some specific links between parameters and observations may provide some empirical limits. The bulge fraction can be reached within a 10% accuracy for all models with the mass ratio of 3 ± 0.5. In the frame of our model for M31, the halo is enriched by stars coming from the progenitors before the interaction, while the thick disk and Giant Stream are further enriched by stars formed before the first passage and fusion. This (see Figure 1 and Section 2) is compliant with pericenter radii of 25 ± 5 kpc. Further modeling and measurements of stellar ages in several fields in the outskirts of M31 would certainly help to provide more accurate values. This would also provide better accuracy for the important epochs of the interaction which would have begun 8.75 ± 0.75 Gyr ago with a fusion time 5.5 ± 0.5 Gyr ago. There is still work needed to properly constraint other orbital parameters such as the inclinations of progenitors relative to the angular momentum direction. We are far from having investigated the whole range of orbital properties which produces a strong tidal tail from stripped material from the satellite, and then a structure similar to the Giant Stream. Table 3 describes all the models presented in this paper, which are part of about 50 tested models. Figure 6 demonstrates that models with differences of progenitor inclinations near 180° show not only a maximum of resonance (thus a very strong tidal tail), but also galactic disks which are probably too warped or distorted. We have mostly investigated models with the main galaxy inclination of 55°–70°, but there could be other solutions for different angles. All the models in Table 3, as well as models in the range described by Table 1 (Column 4) reproduce quite well the observations, including the Giant Stream. A much more detailed analysis is necessary to determine the best model parameters, and this could be done by comparing in detail the constraints provided by the gaseous disk, the thick disk, and the Giant Stream shapes as well as by any other structures discovered in the M31 halo field.

While writing this paper, we have become aware of the discovery of many other structures in the halo of M31 (Mackey et al. 2010), especially after the successful deep imagery of the whole field surrounding M31 by the PANDA team (McConnachie et al. 2010). This magnificent image reveals the presence of loops surrounding M31 which strikingly resemble our model predictions (see Figure 8). In particular, the northern loop found by the PANDA team which extends up to 120 kpc from the galaxy center provides particularly strong constraints on our model. Indeed, the first loop in Figure 8 can accurately reproduce the northern loop (called the NW stream in the incomplete image of Mackey et al. 2010), and other loops are well matched with the eastern arc and Stream C. Figure 11 shows two examples of models that succeed in reproducing most of these features along with the Giant Stream and most other substructures (thin disk, bulge, and thick disk) of M31. In fact, we also verified that the Giant Stream kinematics is still reproduced. Fitting these additional structures allows us to fix the rotation angle around the thin disk axis as well as solving the initial conditions of the system because only the inversion of the initial (x, y, z) coordinate system is able to reproduce the ensemble of features in Figure 11. The size of the loop also depends on the time elapsed after fusion, and the beginning of the interaction is expected to be 8.75 ± 0.35 Gyr within the range of parameters in Table 1.

Zoom In Zoom Out Reset image size

As a whole, M31 and the complex structures in its halo can be reproduced by assuming a single major merger which may have begun 8.7 Gyr ago. The advantages of such a solution for M31 compared to the hypothesis of numerous minor mergers are as follows.

• 1.  
  It explains most of the complexities of M31 by a single event rather than numerous and minor merger events which have to be adapted in a somewhat ad hoc way to each feature discovered in the M31 halo, including the thick disk.
• 2.  
  It naturally explains stellar ages and metallicities of most substructures in the M31 halo as well as why they show so many similarities.
• 3.  
  It overcomes the increasing difficulty of identifying the residuals of satellites assumed to be responsible for the different features in the M31 halo.
• 4.  
  It is consistent with the kinematics of the M31 globular system (e.g., Bekki 2010).

Re-formation of disks is expected after gas-rich major mergers (Barnes 2002; Springel & Hernquist 2005). It is more than just a curiosity since the disk rebuilding scenario has been proposed as a channel for the general formation of spirals (Hammer et al. 2005), including specifically the case of M31 (Hammer et al. 2007). Gas richness of progenitors is certainly high, as it has been estimated for many high-z galaxies. The requirement of a star formation with a low efficiency before the fusion can be reached through any mechanism that prevents star formation in relatively primordial medium. It could be caused either by higher feedback efficiency in primordial medium as assumed in Section 4.2, or because the gas in the progenitors is less concentrated than in present-day spirals as it is in present-day low-surface brightness galaxies, or because cooling is less efficient in a relatively pristine medium, or a combination of all these factors. More generally, the expected increase of the gas metal abundance (expected to be slow before the fusion but very efficient during the fusion) may help to increase the molecular gas fraction, the optical depth of the gas, and the radiation pressure effects, all contributing to a change in the star formation history during the interaction (T. J. Cox 2010, private communication). The above is quite well illustrated a contrario by Bournaud et al. (2010), who assumed that progenitors of high-z major mergers are very strong star-forming, turbulent, gas-rich disks: they naturally find that remnants cannot be disk dominated within such conditions.

There is still a considerable need to improve the modeling of M31, by using a much larger number of particles as well as by providing more kinematical and stellar age constraints in the numerous fields of view surrounding M31. We also need to reproduce with more accuracy the shape of the Giant Stream and of the thick-disk contours. Metal abundance enrichment can also be implemented to more accurately fit the observations. For example, the striking differences between features at different metallicity in the halo found by the PANDA team cannot be retrieved by our modeling in the absence of such an implementation.⁷ We have attempted to verify whether such a trend is predicted by our model. At first order, one may "paint" each star in the progenitors assuming a reasonable metallicity gradient. We adopt the approach by Rocha et al. (2008) assuming that metallicity gradient follows the dust gradient, with a scale length 1.4 times larger than that of the stars. The range of metal abundance (Fe/H) is adopted from 0 to −2.5. The stellar metal abundance of the satellite is likely significantly smaller than that of the main progenitor because their stellar mass ratio (for our 3:1 models, see Table 1) is 5. Little is known about the shape of the metal–stellar mass relationship at z ∼ 1.5, although a serious attempt has been made at z ∼ 0.7 by Rodrigues et al. (2008). Assuming this holds up to z = 1.5, stars in the main progenitor are on average 0.4 dex more metal-rich than those from the satellite, a value smaller than the large scatter of the Rodrigues et al. (2008) relationship. Figure 12 shows the metallicity distribution of "painted" stars in both the Giant Stream (defined as in Section 4.4, from 25 to 80 kpc and Φ_obs from −90° to −140°) and in the northern loop (defined from 55 to 130 kpc and Φ_obs from 0° to 90°) discovered by the PANDA team. Because most of the stars in the northern loop are coming from the satellite, this could explain why this feature is not detected in the metal-rich map of the PANDA team (see Figure 12). Such a prediction may stand unless intermediate-age stars (those formed before the fusion) strongly affect the metal abundance of both the Giant Stream and the northern loop. On the other hand, the same model correctly predicts the (observed) small fraction of such stars in the Giant Stream.

Zoom In Zoom Out Reset image size

Another important improvement in this modeling would be to consider the possible interaction of M33 which has been proposed to be responsible for several features in the SE of the M31 halo by McConnachie et al. (2009). Although the dark matter density profile does not significantly affect our results (results shown in Figures 2–5, 7–9, and 11 are similar when using a Hernquist profile instead of a core profile), we also need to verify whether these are affected by adopting different baryonic mass fractions, especially for lower values than our adopted 20%. Indeed, a considerable change of the dark matter fraction may prevent stars from being stripped by the collision and then the distribution of baryonic matter in the outskirts may be significantly affected.

However, we believe the result of such a model is very encouraging because it fulfills the Occam's razor principle: a single event may explain most—if not all—the properties of M31 and of its outskirts. It also proposes a mechanism for the Giant Stream which is simply fed by stars which are captured into loops orbiting around M31 for elapsed times which may reach a Hubble time. If M31 is actually the result of an ancient major merger, there will be a considerable need to revise our knowledge about our immediate environment and in cosmology,⁸ for the following reasons.

• 1.  
  Merging of gas-rich distant galaxies can easily produce large thin disks by assuming a less efficient star formation before the fusion; this simple scheme supports the rebuilding disk scenario for many giant spiral galaxies, at least for those sharing similar properties (enriched halo stars) as M31.
• 2.  
  Most models show the formation of tidal dwarves, including those from tidal tails formed after the first passage, and some of them could be part of the satellite system of M31; their properties, after comparison with observations, may be used in the future as a constraint to the orbital parameters in further modeling.
• 3.  
  Up to 15% of the material is ejected after the fusion and this may populate the whole Local Group including in the direction of the Milky Way; this led Yang & Hammer (2010) to investigate whether the Large Magellanic Cloud could be related to some ejected material coming from M31.

This work has been supported by the China–France International Associated Laboratory "Origins" and by the National Basic Research Program of China (973 Program), no. 2010CB833000. Part of the simulations have been carried out at the High Performance Computing Center at National Astronomical Observatories, Chinese Academy of Sciences, as well as at the Computing Center at Paris Observatory. We thank Tom Brown and T. J. Cox for their very useful comments during the submission process of this paper. Suggestions and comments from an unknown referee have been very helpful in improving the current version of the article.

Note added in proof: The specific frequency of globular clusters in M31 appears to be three or four times greater than it is in the Galaxy (van den Bergh 2010), i.e., it is expected if a major merger occurred in the first galaxy. The fact that the globular clusters in M31 and in the Galaxy have almost exactly the same half-light radii supports the notion that M31 formed by the merger of two Galaxy-sized ancestors (S. van den Bergh 2010, private communication).