Formation of the Andromeda giant stream: asymmetric structure and disc progenitor

Abstract

We focus on the evidence of a past minor merger discovered in the halo of the Andromeda galaxy (M31). Previous N-body studies have enjoyed moderate success in producing the observed giant stellar stream (GSS) and stellar shells in M31's halo. The observed distribution of stars in the halo of M31 shows an asymmetric surface brightness profile across the GSS; however, the effect of the morphology of the progenitor galaxy on the internal structure of the GSS requires further investigation in theoretical studies. To investigate the physical connection between the characteristic surface brightness in the GSS and the morphology of the progenitor dwarf galaxy, we systematically vary the thickness, rotation velocity and initial inclination of the disc dwarf galaxy in N-body simulations. The formation of the observed structures appears to be dominated by the progenitor's rotation. Besides reproducing the observed GSS and two shells in detail, we predict additional structures for further observations. We predict the detectability of the progenitor's stellar core in the phase-space density distribution, azimuthal metallicity gradient of the western shell-like structure and an additional extended shell in the north-western direction that may constrain the properties of the progenitor galaxy.

galaxies: individual: (M31), galaxies: interactions, galaxies: kinematics and dynamics

1 INTRODUCTION

The Local Group is a natural laboratory for investigating the formation and evolution of galaxies and comparing the observations with theoretical studies. According to the generally accepted cold dark matter model, a snapshot of the Local Group should record a history of the hierarchical structural formation of the Universe. In fact, by studying the spatial, kinematic and metallicity distributions of substructures such as tidal debris and dwarf satellite galaxies in their host halo, we can probe the formation history of galaxies, the density profile of the host galaxy and accretion history of massive black holes (MBHs) associated with satellite galaxies.

Recent deep photometric observations of the halo of the Andromeda galaxy (M31) have discovered a wealth of faint structures, including past and on-going galaxy mergers (Ibata et al. 2001; Irwin et al. 2005; McConnachie et al. 2009; Martin et al. 2013). A giant stellar stream (GSS) and fan-like stellar structures have been discovered in the halo of M31. Its spatial, metallicity and line-of-sight velocity distribution have been observed in detail (Ferguson et al. 2004; Ibata et al. 2004; Kalirai et al. 2006; Gilbert et al. 2009). Wider and deeper surveys of M31's halo region those of Pan-Andromeda Archaeological Survey (PAndAS; McConnachie et al. 2009; Martin et al. 2013) and the Spectroscopic and Photometric Landscape of Andromeda's Stellar Halo survey (SPLASH; Kalirai et al. 2010; Tollerud et al. 2012). In particular, the GSS lies south-east of M31's centre and extends more than 120 kpc along the line of sight (McConnachie et al. 2003; Conn et al. 2016). The fan-like structures spread to the north-east and west side of M31 and with approximate radii of 30 kpc (Fardal et al. 2007). Various works have explored the formation of these structures by colliding two galaxies in N-body simulations, assuming minor (Fardal et al. 2007) or major merger scenario (Hammer et al. 2010). We adopt a past radial interaction model of a dwarf satellite galaxy, which well reproduces almost all of these structures (Fardal et al. 2007, 2013; Mori & Rich 2008; Kirihara, Miki & Mori 2014; Miki et al. 2014; Sadoun, Mohayaee & Colin 2014; Miki, Mori & Rich 2016).

The total mass of the progenitor has been estimated by several approaches. The lower limit of its stellar mass, estimated from the kinematics and luminosity of the GSS, is 10⁸ M_⊙ (Font et al. 2006). Mori & Rich (2008) reported the upper limit of its dynamical mass as 5 × 10⁹ M_⊙. They considered the effect of dynamical friction on the thickness of M31's disc. Fardal et al. (2013) and Miki et al. (2016) also examined the stellar mass of the satellite progenitor assuming a spherical progenitor galaxy. Their best-fitting value was approximately 3–4 × 10⁹ M_⊙.

Fig. 1 shows the stellar mass of the observed dwarf galaxies around M31 versus the projected distance from the centre of M31. The predicted mass range of the progenitor dwarf galaxy is dominated by dwarf ellipticals and dwarf irregulars. On the contrary, all such satellite galaxies in M31 have rotating stellar and/or gas components (McConnachie 2012). Even M32, which is classified as a compact elliptical, has a rotating stellar component with a measured velocity of 55 km s⁻¹ (Bender, Kormendy & Dehnen 1996). M33 is a disc galaxy with no clear bulge. Disc galaxies comprise a large proportion of the less massive galaxies in the local universe (Moffett et al. 2016). Although previous studies have usually assumed a spherical, non-rotating progenitor galaxy, these reasonable conditions motivate us to assume a disc-like galaxy as the GSS progenitor.

Figure 1.

Stellar mass M_* and morphological types of the observed dwarf galaxies surrounding M31 as functions of projected distance from the centre of M31 D_M31(data are taken from McConnachie 2012). Symbols indicate dSphs (open red squares), dIrrs (filled blue circles), dEs (filled magenta triangles), a compact elliptical galaxy (filled cyan square) and an Sb galaxy (inverted filled triangle).

Using star count maps around M31, McConnachie et al. (2003) analysed the surface brightness profile in the direction orthogonal to major axis of the GSS. They obtained an asymmetric surface brightness of the GSS, with sharply decreased star counts at the north-eastern side of the GSS (viewed from the most luminous direction of the GSS). On the other hand, the surface brightness distribution extends widely and smoothly at the western side. This asymmetrical structure of the GSS has never emerged in simulations that assume the infall of a spherical, non-rotating dwarf galaxy (e.g. Fig. 13 d; Gilbert et al. 2007). The characteristic surface brightness profile of the GSS would be an excellent tracer of the morphology of the disrupted progenitor galaxy. To examine the disc merger scenario and identify the progenitor conditions that would lead to complicated evolution, we must systematically scrutinize a large parameter space. Few disc satellite models have generated an asymmetric structure for the GSS, moreover, they have not reproduced the observed shape (Fardal et al. 2008; Sadoun et al. 2014). For example, Fardal et al. (2008) reported two arc-like structures on the eastern side of the GSS (they resemble streams C and D in Ibata et al. 2007). Notably, the sharp edge-like structure at the eastern side of the GSS has not been reproduced.

Owing to its estimated mass, the progenitor galaxy is unlikely to be a nearby dSph, which has mass-to-luminosity ratios of ∼100. Initially, the progenitor is thought to have inhabited a dark matter halo, with a mass-to-luminosity ratio (M/L = 2–5). Similar M/L values have been reported for local satellite dwarf galaxies (see fig. 11 of McConnachie 2012). On the other hand, most previous works have ignored the dark matter halo component of the progenitor galaxy, because this component was rationalized to have been stripped before the first interaction with M31. However, if an inner region of the dark matter halo was gravitationally strongly bound, it would survive the collision with M31. Further it is not understood how the collision would disperse the dark matter halo component through the host halo.

The formation history of galaxies can be inferred from the spatial distribution of heavy elements in their stars. Ibata et al. (2007) obtained the metallicity distribution in the southern area of M31 from colour–magnitude diagrams. They suggested a clear metallicity difference between an eastern high-surface-brightness region (metal rich) and a faint western region (metal poor) of the GSS. Similar trends in the GSS appear in spectroscopic measurements of the metallicity distribution based on the Ca ii triplet absorption lines (Gilbert et al. 2009). In addition, Fardal et al. (2012) measured the metallicity distribution of the western shell along the minor axis of M31's disc. Only recently, Conn et al. (2016) observed the radial metallicity distribution along the GSS. Radial metallicity gradients are commonly observed in dwarf galaxies (Koleva et al. 2009), and are also known in nearby disc galaxies (Magrini et al. 2016). The present-day metallicity gradient of the GSS could conceivably have originated in the progenitor dwarf galaxy (Fardal et al. 2008; Miki et al. 2016).

In surveys of N-body simulations, this paper explores galaxy collisions between M31 and a dwarf satellite galaxy composed of a stellar disc, a stellar bulge and a dark matter halo component. The aim is to reproduce the asymmetric surface brightness of the GSS. In Section 2, we describe the observational data and their treatment, including our simple analysis of the asymmetric surface brightness profile of the GSS. In Section 3, we introduce our modelling of the M31 potential, the N-body satellite progenitor and the numerical model for systematic surveys. The results of numerical simulations and quantitative comparisons with observed data are displayed in Section 4. The metallicity distribution, distribution of the disrupted dark matter halo component of the progenitor, the position of a hypothetical MBH (initially centralized in the progenitor galaxy) and the extended stellar shell at the north-western area of M31 are described in Section 5. We summarize our findings in Section 6.

2 OBSERVED STRUCTURES

2.1 Spatial faint structures around M31

Merger remnants can reveal the properties (mass and morphological type) not only of the host galaxy, but also of the disrupted progenitor galaxy. Some of the faint stellar structures in M31's halo have been well reproduced by an on-going merger of a satellite galaxy (Fardal et al. 2007, 2008, 2013; Mori & Rich 2008; Kirihara et al. 2014; Miki et al. 2014, 2016; Sadoun et al. 2014). Irwin et al. (2005) intensively surveyed the M31 halo with the Isaac Newton Telescope Wide–Field Camera (INT/WFC). Their map covers a 4| $^\circ$ | elliptical region of the semimajor axis with an aspect ratio of 5:3 and an additional ∼10 (deg²) extension towards the south of M31. We analyse their stellar data set in comparisons with our numerical simulations.

The star count maps include non-GSS components such as the halo stars of M31, foreground halo stars of Milky Way and other substructures around M31. To clarify the GSS structure, we simply subtract a constant background as follows. We first calculate the GSS background by removing the area of M31's disc, the eastern and the western shell, the GSS and several stellar substructures from the star count map. Next, we average the star counts per cell, excluding the above non-background area. The stellar background count is 1.18 stars per cell. This count is subtracted from the star counts in each cell. This treatment provides a clearer structure of the GSS.

2.2 Asymmetric structure of the GSS

In their analysis of the star count maps, McConnachie et al. (2003) found an asymmetric surface brightness across the GSS. Specifically, they reported a sharp increase in star counts at the eastern side of the GSS, relative to the western side. To compare our simulated GSS with the observed GSS, we reanalyse the observed data and obtain the azimuthal surface brightness profile of the GSS.

Within the radius R = 2| $_{.}^{\circ}$ |5 (∼30 kpc), the pure GSS component is obscured by superposition of M31's disc and the clumpy stellar structures. Beyond R = 3| $_{.}^{\circ}$ |5, the star counts are insufficient for a proper analysis. Therefore, we analyse the region 2| $_{.}^{\circ}$ |5 < R < 3| $_{.}^{\circ}$ |5 from the centre of M31 over the azimuthal angular range | $30^\circ <\theta _a<100^\circ$ | (east to south). The radial distance is divided into inner and outer areas separated at R = 3| $_{.}^{\circ}$ |0. The analysed area is enclosed by solid lines in Fig. 2(a).

Figure 2.

(a) Background-subtracted stellar count map of the halo of M31 (Irwin et al. 2005). Arcs denote the observed edges of the eastern and western shells and squares locate in the previously observed area of the GSS (McConnachie et al. 2003; Guhathakurta et al. 2006; Fardal et al. 2007). The white ellipse traces the disc of M31. The azimuthal angle θ_a is taken from the eastern direction in M31-centred coordinates. (b) Azimuthal star count distribution of the GSS. The analysis area is a portion of an annulus (⁠| $30^\circ <\theta _a<100^\circ$ |⁠), as shown in panel (a). The dotted and dashed lines profile the inner and outer regions of the annulus, respectively.

Fig. 2(b) shows the azimuthal distribution of star counts in the GSS. Hereafter, we define the brightest azimuthal angle (⁠| $\theta _a\sim 65^\circ$ |⁠) as the GSS axis. The surface brightness gradually decreases at the western side (larger θ_a) of the GSS axis, but increases sharply at the eastern side. Each side of this asymmetric star count profile is fitted by an asymmetrical exponential function proportional to exp| $(-|\text{d} \theta _a|/\delta _{\rm {obs}}^{\pm })$ |⁠, where | $\delta _{\rm {obs}}^{-}$ | and | $\delta _{\rm {obs}}^{+}$ | are the observed widths of the eastern edge and the broad western structure, respectively. The fitted values are summarized in Table 1. The significant width difference between the two sides confirms the strongly asymmetric spatial distribution.

Table 1.

Fitted values of the GSS surface brightness.

Region	GSS axis	Eastern edge	Western broad
	θaxis (⁠| $^\circ$ |⁠)	width | $\delta _{\rm obs}^-$ | (⁠| $^\circ$ |⁠)	width | $\delta _{\rm obs}^+$ | (⁠| $^\circ$ |⁠)
Inner	64.0	2.74	24.1
Outer	65.1	3.34	21.9

3 NUMERICAL MODELS

3.1 M31 potential model

In this work, we assume a fixed gravitational potential model for M31 with a Hernquist bulge (Hernquist 1990), an exponential disc and a dark matter halo with an NFW profile (Navarro, Frenk & White 1996). The scale radius and total mass of the bulge component are 0.61 kpc and 3.24 × 10¹⁰ M_⊙, respectively. The M31 disc is assigned a scaleheight of 0.6 kpc, a radial scalelength of 5.4 kpc and a total mass of 3.66 × 10¹⁰ M_⊙. The inclination and position angle of the M31 disc are | $77^\circ$ | and | $37^\circ$ |⁠, respectively (Geehan et al. 2006). The NFW halo has a scale radius of 7.63 kpc, and a scale density of | $6.17 \times 10^7 \,\mathrm{M}_{\odot }\, \rm {kpc}^{-3}$ |⁠. Under these conditions, the M31 model nicely reproduces the observed surface brightness of the bulge and disc components, the velocity dispersion of the bulge and the rotation curve of the disc (Geehan et al. 2006; Fardal et al. 2007).

Also using full N-body simulations, Mori & Rich (2008) examined the dynamic response of the interaction between M31 and a less massive progenitor (≲5 × 10⁹ M_⊙). They found little change in the gravitational potential of M31. Therefore, M31 can be feasibly treated as a fixed gravitational potential, and the above-mentioned M31 parameters have been adopted in previous studies using the same initial orbital variables of the progenitor (Fardal et al. 2007, 2008; Miki et al. 2016).

3.2 N-body satellite models

To elucidate the origin of the asymmetric GSS surface brightness, we collide M31 with a disc dwarf satellite galaxy, which is a self-consistent N-body disc with stars and dark matter. The initial position and velocity vectors of the progenitor galaxy are taken from Fardal et al. (2007).

An equilibrium model of the disc dwarf galaxy is constructed using in the public code galactics (Kuijken & Dubinski 1995), which generates a self-consistent bulge–disc–halo N-body system. As a tentative disc progenitor model for the THIN model, we adopt the downsized model of the M31 Model A presented in Widrow, Perrett & Suyu (2003), which reduces the M31 mass by a factor of 100. This model treats the density profile of the dark matter halo as a lowered Evans profile (Kuijken & Dubinski 1994), which satisfies an isothermal distribution function with a characteristic radius r_a. The disc component follows an exponential density profile with a scalelength R_d in the radial direction, and an isothermal profile with a scaleheight Z_d in the vertical direction. The bulge component is modelled as a King sphere.

We now summarize the input parameter set of the disc dwarf progenitor models. In galactics, the length, velocity and mass units are 1 kpc, 100 km s⁻¹ and 2.325 × 10⁹ M_⊙, respectively (see Kuijken & Dubinski 1995 and Widrow et al. 2003). The halo parameters input to THIN are the central potential Φ(0) = −1.483, the central velocity dispersion V₀ = 0.952 and r_a = 0.981. The total mass, scalelength and scaleheight of the disc component are M_d = 0.318, R_d = 1.11 and Z_d = 0.13, respectively. The outermost cutoff radius of the disc R_outer is 8.0, with a truncation length of 0.2. The parameters of the bulge component are the central density (ρ₀ = 6.68), velocity dispersion (σ_b = 0.508) and cutoff potential (ψ_cut = −0.835). The cutoff potential controls the extent of the bulge.

As is well known, the nearby dwarf galaxies are relatively thick (Spolaor et al. 2010; Toloba et al. 2011). Accordingly, we construct two additional models (THICK and HOT), in which the scaleheights are Z_d = 0.52 (four times the scaleheight of the disc of THIN) and Z_d = 1.11 (the same length as the scalelength), respectively. The input parameter values are listed in Table 2. The models are seeded with 203 418 particles; 153 752 dark matter particles, 36 756 disc particles and 12 910 bulge particles. When checking the convergence of the numerical resolution, we multiplied this total particle number by five (to 1017 090). Because the galactics code uses a Poisson solver, the total mass in our models is altered by changing the disc thickness. Table 3 lists the ratios of the rotation velocity to the velocity dispersion of the stellar component and the masses of the bulge, disc and dark matter halo in each disc model. The effective radius of the constructed bulge (∼200 pc) is consistent with the relationship between the radius and stellar mass (Gadotti 2009).

Table 2.

Properties of the progenitor models.

Model	THIN	THICK	THICK2	THICK3	THICK4	THICK5	THICK6	THICK7	THICK8	THICK9	HOT
Zd	0.13	0.52	0.52	0.52	0.52	0.52	0.52	0.52	0.52	0.52	1.11
Ψ0	−1.483	−1.483	−1.400	−1.450	−1.500	−1.550	−1.600	−1.650	−1.700	−1.750	−1.483
V0	0.952	0.952	1.000	1.100	1.200	1.300	1.400	1.500	1.600	1.650	0.952

Table 3.

Mass abundances of the progenitor models.

Model	THIN	THICK	HOT
Zd (kpc)	0.13	0.52	1.11
vmax/σ	10.8	2.5	1.3
Bulge (M⊙)	2.9 × 108	3.1 × 108	3.1 × 108
Disc (M⊙)	7.8 × 108	7.3 × 108	6.5 × 108
DM halo (M⊙)	3.2 × 109	3.5 × 109	3.9 × 109

To expand the parameter space of the THICK model, we fix Z_d and vary the rotation curve, which would affect the shape of the GSS. The different rotation velocity of the disc is simply varied by changing the central gravitational potential of the dark matter halo. The resulting models are named THICK2–THICK9 (see Table 2). Fig. 3 shows rotation curves of the various disc models. The THICK disc models are numbered in order of increasing rotation velocity of the disc inner 5 kpc. All of the progenitor models are consistent with the observed baryonic Tully–Fisher relation in the scatter range (McGaugh et al. 2000).

Figure 3.

Rotation curves in various disc models [THICK (lowest disc rotation velocity) to THICK9 (highest disc rotation velocity)].

In the N-body simulations, the gravity is calculated by an original parallel tree code with a tolerance parameter of 0.5. The Plummer softening length is ∼8 pc. These parameters sufficiently resolve the bulge component and reduce numerical two-body relaxation. The orbit is integrated by a second-order leapfrog integrator with a shared timestep of approximately 10 kyr. When testing the convergence of the numerical resolution, we perform an additional high-resolution run using the same code. As mentioned above, the total particle number in the high-resolution run (1017 090) is five times that of the normal resolution. Additionally, we increase the particle number to 16 777 216 (∼16 times that of the high resolution) and conduct a highest resolution run. Besides verifying convergence, the highest resolution run reveals the metallicity distribution in the faint regions such as the broad western structure of the GSS, by which we reduce the Poisson noise in the N-body simulations. In the highest resolution model, we employ the gravitational octree code optimized for graphics processing units (GPU), gothic (Miki & Umemura, in preparation), which implements a hierarchical time step with a second-order Runge–Kutta integrator and the multipole acceptance criterion proposed by Warren & Salmon (1993) and Salmon & Warren (1994). The accuracy control parameter Δ_acc is set to be 2⁻⁸. The resolutions and codes of the normal, high and highest resolution models are summarized in Table 4. Numerical calculations are carried out on the T2K-Tsukuba System, HA-PACS System, COMA System and a workstation at the Center for Computational Sciences, University of Tsukuba, Japan.

Table 4.

Resolution (particle number N) and programming codes at different resolutions.

Resolution	Normal	High	Highest
N	203 418	1017 090	16 777 216
Purpose	Large parameter	Convergence test	High-quality
	survey	for large survey	data analysis
Code	Original tree	Original tree	gothic
	on CPUs	on CPUs	on GPU
Section	Section 4	Section 5.1	Sections 5.2–5.4

3.3 Rotation of satellite progenitor

To reproduce the observed shapes of the GSS and shells, we require an initial inclination of the disc progenitor, which is precluded in spherical progenitor models. Therefore, we systematically vary the initial inclination of the disc and simulate galaxy merger between M31 and the progenitor galaxy. In this survey, we construct a 2D plane in spherical coordinates (ϕ″, θ″) [later altered to | $\phi (\equiv 180^\circ -\phi ^{\prime })$ |⁠, | $\theta (\equiv 180^\circ -\theta ^{\prime })$ |](as shown in Fig. 4). Initially, the pole of the satellite system aligns with the direction of the angular momentum vector of M31's disc. In other words, when | $(\phi ^{\prime },\theta ^{\prime })=(0^\circ ,0^\circ )$ |⁠, the disc planes and the spin axes of the progenitor and M31 discs are identical. The X(Z)-axis corresponds to the minor (major) axis of M31's disc and the arrow points to the north-western (eastern) side of M31's centre. The polar angle θ″ ranges from the pole (+Y) direction to the −X direction (⁠| $0^\circ \leqq \theta ^{\prime }<180^\circ$ |⁠) and the azimuthal angle ϕ″ is measured anticlockwise from the −X-axis on the disc plane in the +Z direction (⁠| $0^\circ \leqq \phi ^{\prime }<360^\circ$ |⁠). The spin axis of the disc dwarf galaxy is inclined first by θ″ and then by ϕ″.

Figure 4.

Initial inclination (ϕ″, θ″) of the spin axis of the progenitor's disc in the specified coordinate system. Grey elliptical disc delineates the M31's disc in the numerical simulation coordinates and the orange elliptical disc is the inclined disc of the progenitor. The green compass shows the north, east and the Local Standard of Rest directions. The Earth locates at the backside of M31's disc. 3D visualization was conducted with the s2plot programming library (Barnes et al. 2006).

Initially, the disc of the progenitor galaxy is inclined by (ϕ, θ), and the (ϕ, θ) plane is comprehensively surveyed with a grid width of | $30^\circ$ |⁠. Next, we carefully examine parameter spaces that reasonably reproduce the GSS. In each disc model (THIN, THICK and HOT), we simulate approximately 350 runs on the (ϕ, θ) plane and identified the appropriate parameter spaces on that plane (e.g. Fig. 11). In the beginning, we explore the case of the models THIN, THICK and HOT, and find an appropriate range of the parameter space on the plane (ϕ, θ) (e.g. Fig. 11). We then perform detailed simulations in the appropriate parameter spaces (⁠| $-50^\circ \le \phi \le 60^\circ$ | and | $10^\circ \le \theta \le 105^\circ$ |⁠). All models are iterated through 74 runs except the THICK7 model (176 runs).

4 NUMERICAL RESULTS

4.1 Representative models

In this section, we simulate galaxy collisions between M31 and a disc dwarf satellite galaxy. First we demonstrate a successful model that well reproduces the GSS axis, width of the sharp eastern edge and width of the broad western structure of the GSS. Fig. 5 shows the whole spatial distribution of the satellite galaxy remnant and the normalized stellar number count, as function of the azimuthal angle for one of the successful parameters [model THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |]. To ensure a high-quality analysis, we here use the data of the highest resolution runs (see Table 4). The shapes and sizes of the north-eastern and western shells are also well reproduced. The westernmost side of the GSS exhibits an arm-like structure rather than the observed smooth structure. However, the true structure is poorly understood because the GSS features are faint in that region and substructures have been observed there. The northern side of the eastern shell presents a dense region at the approximate position of the northern spur reported by Ferguson et al. (2002). The faintest structures are an extended shell structure outside the western shell (Section 5.4) and a spillover at the tip of the GSS (⁠| $\eta \sim -5^\circ$ | and | $\xi \sim 4^\circ$ |⁠). Ibata et al. (2007) reported a similar structure called Stream B, which is mainly composed of metal-rich stars. Whereas the observed structure is almost perpendicular to the GSS, the simulated one has a smaller interior angle.

Figure 5.

(a) Surface mass-density distribution of the disrupted progenitor, under the parameters yielding one of the most successful results [model THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |]. The inclined elliptical line describes the shape of M31's disc. Square symbols indicate the observed fields of the GSS (McConnachie et al. 2003; Guhathakurta et al. 2006). Circles are the edge positions of the eastern and western shells analysed by Fardal et al. (2007). (b) Normalized stellar count in the GSS as a function of azimuthal angle. Blue solid and black lines present the N-body simulation results and the observed profile of the GSS, respectively.

Fig. 6 shows the time evolution of the disrupted dwarf galaxy on the sky coordinate up to current epoch. The progenitor galaxy, initially located in the north-western area of M31 (Fig. 6i), collided almost head-on with M31. From its initial condition, the progenitor disc reached the first pericentric passage of the orbital motion in approximately one dynamical time of the progenitor's disc. The progenitor is then disrupted and a component spreads into the south-eastern area with the growth of the GSS (Fig. 6j). Simultaneously, part of the debris falls and enters the western side of M31's centre. After the second pericentric passage, the debris spreads into a wide fan called the eastern shell (Fig. 6k). Immediately, some of the debris moves to the western area, forming a similar shell called the western shell (Fig. 6 l). In each run, the current epoch is defined as the snapshot, in which the simulated edge positions of the eastern and western shells best match the observed positions.

Figure 6.

Evolution of the surface mass-density distribution of the disrupted dwarf galaxy in sky coordinates. Top to bottom panels show the evolution of the bulge component (a–d), disc component (e–h), bulge and disc components (i–l) and dark matter halo component (m–p). Left-hand to right-hand panels present the mass-density distributions at 0.0, 0.3, 0.6 and 0.9 Gyr (current epoch) after the start of the simulation run. Symbols and lines in each panel are those of Fig. 5(a).

In Fig. 6(l), the dense region of the simulated GSS lies along the observed fields (indicated by open squares; McConnachie et al. 2003; Guhathakurta et al. 2006). In the simulation results, the boundaries of the shell structures at the north-eastern and western areas reach the observed edges of the shells analysed by Fardal et al. (2007, indicated by black circles). Most of the bulge component resides in the eastern shell (Fig. 6d); almost none is found in the GSS area. This occurs because the bulge stars are strongly bound by the gravitational potential of the progenitor's bulge and therefore survive the tidal disruption of M31's gravitational field at the first pericentric passage. The simulated progenitor also exhibits a spherical dark matter halo component, which evolves as shown in the bottom panels of Fig. 6. The spatial distribution spread over a quite broad region around M31, especially at the eastern side of the GSS. Interestingly, the sharp edge of the dark matter distribution in the western M31 appears to locate at the stellar component.

Fig. 7 shows the mass fraction of the dark matter component originally surrounding the progenitor dwarf, relative to M31. The assumed mass distribution of the dark matter halo inhabited with M31 is described in Section 3.1. This figure is constructed from the snapshot of Fig. 6(p). The highest mass fraction is approximately 0.1. Relative to the smooth dark matter halo of M31, the progenitor's dark matter is significantly enhanced at the edge of the western shell. Interestingly, this region of mass-density enhancement coincides with the stellar structure at the western shell. However, this result indicates that the mass fraction is too small to detect even with thirty-metre telescope, which exploits the weak lensing of the background halo stars.

Figure 7.

Mass ratio of the dark matter halo initially inhabited with the progenitor galaxy, relative to the dark matter halo of M31. The inclination of the disc model is | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ | (model THICK7).

Fig. 8 displays the 3D distribution of the debris in the spherically symmetric and axisymmetric progenitor models. The merger with a spherical galaxy is simulated under the same mass-resolution (11 520 000 particles) as the baryonic component of the disc model in the highest resolution run. The spherical dwarf satellite galaxy has a Plummer equilibrium distribution with a total mass of 2.2 × 10⁹ M_⊙ and a scale radius of 1.03 kpc. The 3D distribution of the simulated GSS is consistent with the latest observations by Conn et al. (2016). The simulated GSS is slightly shorter than the observed GSS, possibly because we imposed an artificial radial cutoff in the initial progenitor's disc for simplicity. Dynamically cold components (e.g. fine structures in the inner regions of the eastern and western shells) appear in the disc merger scenario (Fig. 8). The most important difference is the GSS width; in particular, a dynamically cold component on the western side of the GSS. The southern and western spread of the GSS is much broader in the disc merger than in the spherical progenitor merger in the 3D view.

Figure 8.

3D distribution of the colliding satellite galaxy in the disc model (both disc and bulge stars are plotted) and the spherical progenitor model at the best-fitting epoch. Symbols and lines in each graphic are those of Fig. 5(a). Observed distances in each region are indicated as follows: open red circles with error bars (best parameters in Conn et al. 2016), filled black triangles (most likely parameter values in Conn et al. 2016) and cyan crosses (McConnachie et al. 2003). Left-hand panels (a–c): results of the successful disc model (progenitor model using the same parameters as Fig. 6 l). Right-hand panels (d–f): results of the spherical symmetric Plummer model. Panels (a) and (d): view on the sky coordinates. Panels (b) and (e): View on the line-of-sight depth d_M31 (kpc) versus η (deg) plane. Panels (c) and (f): view on the line-of-sight depth d_M31 (kpc) and ξ (deg) plane. White arrows show the line-of-sight direction from Earth.

Fig. 9 shows the 3D distributions of the bulge and the dark matter halo components at the best-fitting epoch. In panels (b) and (c), the progenitor's bulge is elongated along the line-of-sight direction by the tidal force of M31's potential. Several shell-like structures appear in the 3D view of the dark matter halo component of the disrupted progenitor, indicating that part of the dark matter halo component crosses the central region of M31 several times.

Figure 9.

3D distribution of the bulge (left-hand panels) and dark matter (right-hand panels) components of the disrupted dwarf galaxy at the best-fitting epoch. Symbols and lines in each panel are those of Fig. 5(a). Viewing angle is the same as in Fig. 8. The upper right graphic in panel (a) enlarges the 1| $_{.}^{\circ}$ |5×1| $_{.}^{\circ}$ |5 region outlined in green. The phase-space distribution of the bulge component is analysed in the white-outlined region (see Section 5.3). The red circle indicates the position and size of a recently discovered density enhancement on M31's disc (Davidge 2012).

We now present some typical results of clockwise- and anticlockwise-rotating disc models in the sky coordinate system. Figs 10(a) and (b) show the stellar mass-density distributions of almost clockwise (⁠| $(\phi ,\theta )=(-90^\circ ,100^\circ )$ |⁠) and anticlockwise (⁠| $(\phi ,\theta )=(90^\circ ,75^\circ )$ |⁠) rotating disc models, respectively. The shape of the debris largely differs from that of the galaxy–spherical dwarf merger. The progenitor–M31 collision was almost head-on, and the progenitor centre passed barely east of M31's centre. The shortest distance between the progenitor and M31 centres was 1 kpc at the first pericentric passage. In addition, the progenitor's disc is more than 1 kpc across (Section 3.2). As the details depend on the inclination of the progenitor disc, the main component of the progenitor passes to the east of M31's centre, but a portion enters the western area of M31's centre.

Figure 10.

Projected stellar mass-density distributions of (a) an initially clockwise-rotating disc model [model THICK with | $(\phi ,\theta )=(-90^\circ ,100^\circ )$ |] in the sky at the best-fitting epoch and (b) an initially anticlockwise-rotating disc model [model THICK with | $(\phi ,\theta )=(90^\circ ,75^\circ )$ |]. Symbols and lines in each panel are those of Fig. 5(a).

Clearly, there is no GSS-like component in Fig. 10(a). Instead, we observe a curious structure at the south-eastern area of M31. This occurs because the stellar component passing to the west of M31's centre at the first pericentric passage has a clockwise-rotating component, so is ejected eastward. On the other hand, in Fig. 10(b), the stellar component passing just east of M31's centre has an anticlockwise-rotating component, so the GSS is slight on the eastern side but spreads broadly on the western side. However, the direction of the simulated GSS offsets relative to that of the observed GSS.

To summarize the above findings, the direction and shape of the GSS strongly depend on the rotation of the progenitor galaxy. To evaluate how the initial inclination of the progenitor affects the GSS, we analyse the reproducibility of the GSS axis, the width of the eastern edge and the width of the broad western extent in the following subsections.

4.2 GSS axis

The direction of the GSS axis explicitly informs the reproducibility of the GSS. As stated in Section 4.1, the azimuthal angle of the density peak in the GSS largely depends on the initial inclination of the progenitor's disc. To examine the effect of this initial parameter, we require a complete sweep of the large parameter space.

To quantitatively compare the observed and simulated GSS axes θ_axis, we conduct a χ² analysis. Mimicking the observed data analysis, we obtain the simulated GSS axis by an asymmetric exponential fitting of the azimuthal distribution of the GSS (as described in Section 2.2). The top panels in Fig. 11 show the χ² maps of the simulated GSS axis at the best-fitting epoch on the (ϕ, θ) plane in each model. Here we assume an observed uncertainty | $\sigma _{\rm obs}^{\rm peak}$ | of | $\delta _{\rm obs}^-$ |⁠, the observed width of the eastern edge (see Table 1). The bluer region well matches the observed GSS direction and the (ϕ, θ) parameter space that reproduces the observed GSS axis is tightly constrained. The thick black line in each panel describes the 1σ confidence interval of Δχ² (=1). The minimum χ² value in the THIN, THICK and HOT disc models are 0.01, 0.28 and 0.44, respectively.

Figure 11.

Top panels: comparisons between the azimuthal angles of the density peaks in the observed GSS and the simulation runs in each disc model on the plane (ϕ, θ) of the initial inclination of the progenitor's disc (χ² analysis). Middle panels: comparisons between the observed and simulated eastern edge widths. Coloured squares show the χ² values of the simulated parameters. Thick solid, thin solid and dashed contour lines indicate the 1σ, 2σ and 3σ confidence intervals of the Δχ², respectively. Thick red contours describe the 1σ confidence intervals of the azimuthal angles of the density peaks in the GSS (top panels). Bottom panels: the initial spin axis of the progenitor's disc on the observed frame is related to the eastern edge width of the GSS. The colour bar shows the inner product of the normalized line-of-sight vector into M31's centre and the normalized initial spin vector of the progenitor's disc. Magenta (green) areas indicate clockwise (anticlockwise) rotation of the progenitor's disc. Within the yellow area, the disc is viewed almost edge-on from the Earth.

The well-fitting area apparently shifts to larger θ with increasing scaleheight of the disc. The reason for this trend is discussed in Section 4.5. Here the χ² value in the Plummer model is 2.5, reasonably suitable for detecting the GSS axis.

4.3 Eastern edge of the GSS

The star counts sharply decrease from the GSS axis in the eastward direction, relative to the western direction. We analyse the eastern edge similarly to the GSS axis.

To compare the observed and simulated widths of the eastern edge, we conduct a χ² analysis assuming an observational uncertainty | $\sigma _{\rm obs}^{\rm east}\equiv \delta _{\rm {obs}}^-$ | (the observed width of the eastern edge stated in Table 1). Mimicking the observed data analysis, we obtain the width of the simulated GSS by an asymmetric exponential fitting. The middle panels in Fig. 11 show the simulated widths of the GSS eastern edge at the best-fitting epoch on the (ϕ, θ) plane in each disc model. The minimum χ² values in the THIN, THICK and HOT disc models are 0.01, 0.04 and 0.18, respectively. Note that in the THIN and THICK models, a specific parameter range replicates both the observed GSS axis and the eastern edge of the GSS.

At present, how the surface brightness profile of the GSS becomes asymmetric is unclear. Here, we examine the mechanism that forms the eastern edge of the GSS on the plane (ϕ, θ). The bluer region well fits the observed width of the GSS eastern edge and a clear boundary of χ² values separates the well-fitting parameters from the remaining (ϕ, θ) parameter space. The colour maps in the bottom panels of Fig. 11 show the directions of the initial inclination of the progenitor's disc, viewed from the Earth. These maps are overlaid on the χ² maps of the eastern edge width of the GSS. The colour bar indicates the inner product of the normalized line-of-sight vector into M31's centre and the unit vector of the progenitor's disc spin. This figure describes the behaviour of the rotating disc in the sky, providing intuitive information on the successful conditions. The boundaries in the middle panels of Fig. 11 resemble the curve of the edge-on region in the sky coordinates. Viewed from the Earth, the spin vector of the progenitor's disc switches between positive and negative as it crosses the curve. In other words, the curve is the switching line of the clockwise and anticlockwise rotation of the progenitor's disc in the sky.

Although the clockwise -rotating disc models nicely replicate the GSS axis, they fail to reproduce the sharp eastern edge. The successful parameter space resides in the plane [(ϕ, θ) with | $-45^\circ <\phi <30^\circ$ | and | $30^\circ <\theta <80^\circ$ | (THIN model) or | $-20^\circ <\phi <20^\circ$ | and | $40^\circ <\theta <70^\circ$ | (THICK model)]. The Plummer model fails to construct the edge and the minimum χ² value is 9.9. This confirms that spherical progenitor models cannot reproduce the sharp eastern edge.

4.4 Broad western structure of the GSS

In contrast to the eastern side, the star counts at the western side of the GSS gradually decrease from the GSS axis to the westward direction. We investigate the internal structure of this region and the reproducibility of the broad western structure, comparing the simulated results with the observed data. For comparison, we set the half width of the broad western structure | $\delta _{\rm {obs}}^+/2$ | (see Table 1) as the uncertainty | $\sigma _{\rm obs}^{\rm {west}}$ | in the χ² analysis. This large uncertainty | $\sigma _{\rm obs}^{\rm {west}}$ | reflects the faint and noisy nature of the wider western region of the GSS. Consequently, the observed distribution fitting might overestimate the true width of the broad western structure. Nevertheless, the western side of the GSS is significantly wider than the eastern edge in reality.

Fig. 12 shows the χ² map of the western side of the GSS at the best-fitting epoch in the (ϕ, θ) plane. The bluer region well reproduces the observed broad western width of the GSS. The minimum χ² value in the THIN, THICK and HOT models are 0.74, 0.97 and 0.21, respectively. Again, we estimate the width of the simulated GSS by an asymmetric exponential fitting. This analysis does not apparently limit the parameter space. In general, the thicker the disc model, the better the reproducibility of the observed western width. The best-fitting parameters in the THICK and HOT models reproduce all three stellar distributions of the GSS: the GSS axis, the eastern edge and the broad western structure.

Figure 12.

Same as Fig. 11, but comparing the observed and simulated western widths of the GSS.

Fig. 13 plots the normalized number count as a function of the azimuthal angle for | $(\phi ,\theta )=(5^\circ ,40^\circ )$ | in various disc and spherical progenitor models. The selected parameters in Fig. 13 yield successful results in the THICK model. In the THIN model, the azimuthal angle of the density peak in the GSS matches the observed one, but the stream is dynamically too cold and is too narrow. In contrast, the THICK model (Fig. 13b) well reproduces the observed axis and eastern edge of the GSS. The western side of the simulated GSS is wider than the eastern edge. At this inclination, the profile of the HOT model resembles that of the spherical model (Fig. 13d), with an angular shift of the GSS axis. The western and eastern sides of the simulated GSS are almost symmetric about the GSS axis. In the Plummer model, the χ² value of the western width is 3.2.

Figure 13.

Azimuthal angle distributions (red lines) of the GSS simulated by the (a) THIN, (b) THICK, (c) HOT and (d) Plummer models. The inclination of the disc models is fixed at | $(\phi ,\theta )=(5^\circ ,40^\circ )$ |⁠. The observed GSS distribution (black dashed lines) is also plotted in each panel.

4.5 Effect of rotation velocity

As stated in Section 4.3, the width of the eastern edge can be explained by the anticlockwise rotation of the progenitor's disc in the sky. As displayed in Fig. 11, the parameters that well fitted the GSS axis are substantially limited on the (ϕ, θ) plane. When around ϕ = 0, the parameter space that well reproduces the edge width favours smaller θ. The THIN model reproduces the GSS axis at smaller θ than the larger scaleheight models. This tendency might be attributable to the varying thicknesses of the disc models. In fact, as mentioned in Section 3.2, changing the thickness of the disc model alters the total mass of the progenitor model, and hence the rotation velocity of the disc model. As varying the rotation velocity of the progenitor's disc would shift the simulated GSS axis, we expand the parameter space of the THICK model (maintaining constant Z_d) to vary the rotation curves.

Fig. 14 presents the results of the additional simulations in and around the acceptable parameter range of the THIN and THICK models. The number assigned to each disc model indicates its rotation speed (lowest for THICK, highest for THICK9; see Table 2). Higher rotational velocity reduces the θ that reproduces the GSS axis. As the rotation speed of the disc increased, the fitted parameter space converges towards the THIN's parameter space (Fig. 11 top-left panel).

Figure 14.

χ² analysis of the azimuthal angle of the GSS density peak in the thick disc models. The symbols are explained in Fig. 11.

Here, we briefly summarize our findings. A disc progenitor with higher rotation velocity shifts the GSS axis θ to smaller values. The eastern edge of the GSS forms only under anticlockwise rotation of the progenitor dwarf, regardless of its disc thickness. In the extremely thin disc model, the debris is too dynamically cold to reproduce the observed spatial extent of the GSS. On the other hand, the HOT model could not form the observed sharp eastern edge. The broad structure at the western side of the GSS likely originates from an anticlockwise rotating component passing just east of M31's centre at the first pericentric passage.

4.6 Line-of-sight velocity distribution

Additionally to the photometric survey of stars, spectroscopic measurements of the line-of-sight velocity distribution in the GSS have also been carried out (Ferguson et al. 2004; Kalirai et al. 2006; Koch et al. 2008; Gilbert et al. 2009). The observed area almost follows the extending direction of the GSS (Ibata et al. 2004). Fig. 15 shows the line-of-sight velocity distribution of the simulated GSS at the best-fitting epoch. The simulated data are those of the highest resolution run described in Section 4.1 [model THICK7; | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |] and the simulation area is | $30^\circ <\theta _a<100^\circ$ | (see Fig. 2a). The distance and heliocentric systemic line-of-sight velocity of M31 are assumed as 780 kpc and −300 km s⁻¹, respectively (Font et al. 2006). The simulated line-of-sight velocity distributions are consistent with that of the observed data.

Figure 15.

Mass density of the line-of-sight velocity distribution of the simulated GSS [model THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |]. Cyan symbols are observed data in each field (Ferguson et al. 2004; Kalirai et al. 2006; Gilbert et al. 2009). Each bar on the symbols indicates the line-of-sight velocity distribution of the stars at that distance.

Gilbert et al. (2009) reported an additional cold component (R ≲ 20 kpc) near the eastern edge of the GSS. This component is absent on the sky coordinates in our result (Fig. 8a), but a similar structure overlaps on the eastern shell [see the 3D map ((ξ, d_M31) = (0| $_{.}^{\circ}$ |5, 5kpc)) in Fig. 8c]. Miki et al. (2016) reported a similar component (third shell) in many parameter sets.

5 DISCUSSIONS

5.1 Availability of the model and assumptions

To test the convergence of the numerical resolution, we perform a high-resolution run using the parameters that well reproduce the observed structures [THICK7; | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |]. The high-resolution run simulates 1017 090 particles (five times the particle number in normal-resolution runs). Fig. 16 shows the convergence results. Panel (c) presents the results of the highest resolution run, which the number of particles is about 16 times higher than in the high-resolution run. All three runs yield similar global structures, but the distribution is somewhat noisy in Fig. 16(a). The azimuthal angle distributions of the GSS are consistent in the normal- and higher resolution runs (Fig. 16d). Therefore, we consider the resolution of our systematic surveys is sufficient.

Figure 16.

Convergence test of our simulations with a successful parameter set [model THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |]. Spatial distributions of the disrupted progenitor at the best-fitting epoch with different resolutions: (a) normal-resolution simulation in parameter surveys, (b) high-resolution simulation, which uses five times more particles than the normal-resolution runs, (c) highest resolution simulation, which is implemented in different code (gothic). Symbols and lines in panels (a), (b) and (c) are those of Fig. 5(a). (d) Azimuthal angle distribution of the GSS. The dotted magenta, solid red and solid blue lines are generated at normal, high and highest resolutions, respectively. The black dashed line denotes the observed data (Irwin et al. 2005).

Miki et al. (2014) systematically evaluated the infalling orbit of a spherical progenitor galaxy. They limited the possible orbital parameters within a narrow range including the orbit proposed by Fardal et al. (2007). The tight constraint originates from the strength of the tidal force exerted by M31's bulge and the passage duration of M31's central region. To form the GSS and both shells, the progenitor must be almost entirely disrupted. Therefore, assuming that the feasible orbital parameters are independent of the progenitor's morphology, we here adopt the orbit found by Fardal et al. (2007). Our hypothesis will be tested in a future study (i.e. a systematic orbit survey of the disc progenitor).

Kirihara et al. (2014) examined the outer density distribution of the dark matter halo of M31. They highlighted the necessity of reproducing the observed surface brightness ratios among the GSS and both shells. They suggested that the varying gravitational potential of M31 changes the evolutional time-scale of the merger remnants and forms appropriate structures. However, in a spherical progenitor merger, their best-fitting parameter did not replicate the characteristic asymmetric structure of the GSS. In addition, their results might depend on the morphology of the progenitor. For this reason, we first assume the generally adopted conditions in our M31 model. On the other hand, our successful disc model does not explain the observed surface density ratio. Therefore, the mass-density profile of M31's dark matter halo, the orbital initial conditions and the detailed morphology of the progenitor are interesting future investigations.

5.2 Metallicity distribution

The GSS exhibits non-uniform metallicity in the perpendicular direction to the GSS axis (Ibata et al. 2007; Gilbert et al. 2009). The observed mean metallicities in the ‘core’ region (high-brightness region; outlined by the green dashed line in Fig. 17) and the ‘cocoon’ region (western envelope; outlined by the magenta dashed line in Fig. 17) are 〈[Fe/H]〉  = −0.54 and 〈[Fe/H]〉  = −0.71, respectively. On the other hand, nearby dwarf galaxies exhibit radial metallicity gradients of −0.6 < Δ[Fe/H] <0.2 (Koleva et al. 2009). The radial metallicity is computed by

\begin{equation} \Delta {\rm [Fe/H]}\ \equiv \frac{{\rm d}\, \rm [Fe/H] ({\it r})}{{\rm d}\, \rm {log}({\it R}/{\it R}_{\rm {d}})}, \end{equation}

(1)

where R_d is the scalelength of the progenitor's disc. The spatial metallicity distribution of the merger remnants could reveal the initial metallicity gradient of the progenitor galaxy. In fact, Fardal et al. (2008) showed that disc infalling models generate the differences of metallicity in the GSS, although the initial radial gradient of the progenitor is quite high [Δ[Fe/H]  ∼ −1.0, read by eye from fig. 2(b) in Fardal et al. 2008]. Metallicity gradients in the GSS appear even in models of spherically symmetric progenitor galaxy (Miki et al. 2016). However, the observed large differences in mean metallicity do not easily develop in these models.

Figure 17.

(a) Selected analysis regions. The background is the grey-scale mass-density distribution of the disrupted progenitor galaxy [model THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |]. (b) Metallicity distribution with a metallicity gradient of Δ[Fe/H]  = −0.5 and the mean metallicity of 〈[Fe/H]〉  = −0.57. (c) Metallicity distribution filtered by the detection limit of INT/WFC. Green and magenta dashed lines, respectively, denote the ‘core’ and ‘cocoon’ regions (Ibata et al. 2007). Other symbols and lines in each panel are indicated in Fig. 5(a).

To estimate the metallicity distribution in faint regions, such as the broad western structure of the GSS, we need to reduce the Poisson noise in the N-body simulations. For this purpose, we seed the progenitor galaxy with over 16 million particles in the highest resolution simulation using gothic (see Table 4). Other parameters for the progenitor model are those of the high-resolution run.

As shown in the top panels of Fig. 6, most of the disrupted bulge components appear on the M31 disc. Therefore, we set the metallicity gradient only to the disc component of the progenitor galaxy. Initially, we assume Δ[Fe/H]  = −0.5, the observed metallicity gradient of Koleva et al. (2009). We also set the mean metallicity to 〈[Fe/H]〉  = −0.57, consistent with the mass–metallicity relation of nearby dwarf galaxies (Dekel & Woo 2003).

Fig. 17(b) shows the metallicity distribution in the disrupted disc of the progenitor galaxy. To remove extremely faint structures from this figure, we set a simple detection limit and generate Fig. 17(c). The total absolute magnitude of the GSS is set to M_V = −14 (Ibata et al. 2001) and the apparent magnitude limit in the V band is 24.5 (detection limit of INT/WFC). The initial radial metallicity gradient yields large difference in the metallicity distribution in the azimuthal direction of the GSS. Fig. 17(c) suggests that in the east–west direction, the stellar population of the GSS originates from the centre of the initial satellite progenitor's disc and proceeds towards the outside. Interestingly, similar azimuthal differences of metallicity also occurred in the western shell (see Fig. 17c).

Ibata et al. (2007) and Gilbert et al. (2009) observed the mean metallicity in the cocoon region, which is far from the ‘core’ of the GSS and contains few simulated particles. Therefore, we cannot directly compare the simulation results with observations in this region. Instead, we analyse the azimuthal metallicity distribution near the radius of the observed data (3| $_{.}^{\circ}$ |5 < R < 4| $_{.}^{\circ}$ |5 from M31's centre). Fig. 18 shows the azimuthal distribution of the mean metallicity of the GSS, where the mean metallicity of the whole progenitor's disc is −0.5, and the initial metallicity gradient is varied from −0.5 to −0.2. The GSS axis is | $\theta _a\sim 65^\circ$ | (see Table 1). The mean metallicity around this axis is relatively high and almost equals the mean metallicity of the whole progenitor's disc. On the other hand, the mean metallicity in the GSS envelope (⁠| $\theta _a\gtrsim 80^\circ$ |⁠) is relatively low. We obtain strong metallicity differences from the core to the envelope regions of the GSS.

Figure 18.

Mean azimuthal metallicity distribution in the GSS. The inclination of the disc model is THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |⁠. The mean metallicity of the progenitor's disc 〈[Fe/H]〉 is set to −0.5. The initial metallicity gradient of each line is −0.2 (red line with circles), −0.3 (blue line with triangles), −0.4 (magenta line with squares) and −0.5 (black line with diamonds).

Fardal et al. (2012) observed the mean metallicity in the western shell along the minor axis of the M31 disc and obtained 〈[Fe/H]〉  = −0.7. Here, we roughly fit our simulated mean metallicity to the observed 〈[Fe/H]〉 values under two metallicity gradient conditions. For the smaller Δ[Fe/H]  = −0.3 and larger Δ[Fe/H]  = −0.5 metallicity gradients, the best-fitting 〈[Fe/H]〉 values are −0.62 and −0.57, respectively.

We also analyse the azimuthal metallicity distributions at similar radii of the observed data. Fig. 19 plots the azimuthal and radial metallicity distributions of the GSS and the western shell for (Δ[Fe/H], 〈[Fe/H]〉) = (−0.5, −0.57) and the case of (Δ[Fe/H], 〈[Fe/H]〉) = (−0.3, −0.62). Fig. 17(a) shows the analysed areas of the mean metallicity distributions in the azimuthal and radial directions (outlined by red solid and blue dashed lines, respectively) of the GSS and western shell. The azimuthal mean metallicity distributions of the GSS are plotted in Fig. 19(a). The clearest metallicity differences appear at 3| $_{.}^{\circ}$ |5 < R < 4| $_{.}^{\circ}$ |5. In the higher metallicity gradient model, the difference of metallicities in the core region (⁠| $\theta _a=65^\circ$ |⁠) and envelope region (⁠| $\theta _a=85^\circ$ |⁠) differs by approximately 0.25 dex. Conversely, as shown in Fig. 19(b), the radial metallicity differences in the GSS are small in narrow azimuthal ranges (⁠| $60^\circ <\theta _a<70^\circ$ | and | $70^\circ <\theta _a<80^\circ$ |⁠). Radial metallicity differences in the GSS were only recently reported by Conn et al. (2016). Although their data are azimuthally averaged, our results are qualitatively consistent with theirs; namely, the metallicity differences are higher near the root than in the tail of the GSS. In Fig. 19(b), we show the azimuthal metallicity differences inside and outside of the western shell. The metallicity differs by approximately 0.2 from south to north inside the shell, but scarcely differs outside the shell. Along the minor axis of M31's disc, the mean metallicities are similar inside and outside the western shell. Fardal et al. (2012) measured only the directional mean metallicities, and their data set stacks the metallicities of stars inside and outside the western shell. Fig. 19(d) shows the radial metallicity differences in the western shell. The mean metallicity drastically decreases at the edge of the western shell (⁠| $155^\circ <\theta _a<175^\circ$ |⁠), suggesting largely inhomogeneous metallicity distribution in the western shell.

Figure 19.

Mean metallicity distributions in (a) azimuthal direction of the GSS, (b) radial direction of the GSS, (c) azimuthal direction of the western shell and (d) radial direction of the western shell. The inclination of the disc model is THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |⁠. Solid and dotted lines are the results of (Δ[Fe/H], 〈[Fe/H]〉) = (−0.5, −0.57) and (Δ[Fe/H], 〈[Fe/H]〉) = (−0.3, −0.62), respectively. The vertical black solid line and the green square in panel (c) indicate the minor axis of M31's disc and the observed mean metallicity along the minor axis (Fardal et al. 2012), respectively.

5.3 Progenitor's bulge and central MBH

We now discuss the current position of an MBH initially centred at the progenitor galaxy of the GSS. According to the hierarchical structure formation of the Universe, MBHs centralized in dwarf galaxies should be wandering in the halo of their host galaxy. The assumed spherical component of the GSS progenitor galaxy has an approximate stellar mass of 10^8–9 M_⊙. The observed mass of a central MBH correlates with the mass and velocity dispersion of the spherical component (Magorrian et al. 1998; Gültekin et al. 2009). The velocity dispersion of the bulge σ_bulge is ∼50 km s⁻¹ and the MBH mass is simply estimated as 4 × 10⁵ M_⊙ assuming the M–σ relation (Gültekin et al. 2009). Baldassare et al. (2015) presented a small central MBH, which has a mass of ∼5 × 10⁵ M_⊙, lies on the correlation. The bulge component of our disc progenitors is somewhat less massive than in the spherical model assumed in Miki et al. (2014). Although not demonstrated in the present simulation, it is important to note that the mass of the bulge component can change by several factors when varying the bulge–disc ratio and the mass-to-luminosity ratio.

Our simulations can track the current position of the putative MBH initially centred in the progenitor. In fact, Miki et al. (2014) predicted the current position of the progenitor MBH in the halo of M31, varying the orbits of the progenitor galaxy of the GSS. The radiation spectrum of the gas surrounding the MBH was estimated from the advection-dominated accretion flow model (Kawaguchi et al. 2014), which reasonably describes this phenomenon. According to their results, the MBH might be observed with currently operating radio band telescopes such as JVLA and ALMA. They assumed a spherical progenitor galaxy of the GSS. We emphasize that the position of the surviving bulge component (Fig. 6d) approximates the most suitable orbit (ID 1) of Miki et al. (2014) in the sky coordinates. Despite the different morphology of our progenitor, its bulge core component at the best-fitting epoch approximates the current position of the MBH predicted by Miki et al. (2014). The bulge component and MBH undergo similar orbital motions under the gravitational potential of M31, which essentially controls both orbits. In addition, the progenitor is currently passing through the apocentre of its orbital motion with a slow drift velocity, implying that the current bulge position is relatively long term. Also, the best-fitting epoch of the simulation runs, as defined by the edges of the eastern and western shells, is less variable than when a spherical progenitor is assumed. Therefore, our results are consistent with the current MBH position predicted by Miki et al. (2014).

One might expect that the progenitor's bulge can trace the MBH and appear in surface brightness data and/or phase-space mass-density distributions. For instance, Davidge (2012) located a north-eastern clump at (ξ, η) = (0.24, 0.20) on M31's disc with an effective radius of ∼600 pc. The position and size of this clump are indicated by the red circle in Fig. 9(a). Davidge (2012) estimated the stellar mass of the clump as 3 × 10⁸ M_⊙, consistent with the mass of the progenitor's bulge in our models (∼3 × 10⁸ M_⊙).

We analyse the phase-space distribution of the progenitor's bulge in the observed frame. We construct an M31 model with N-body particles in the MAny-component Galactic Initial-conditions generator (Miki & Umemura in preparation), adopting the physical quantities of our present fixed potential model. Fig. 20 shows the phase-space mass-density distributions of the disrupted progenitor galaxy and M31 stars. This figure is constructed from the same data as Fig. 16(c). The analysed region is outlined in white in Fig. 9(a). The clumpy region at R ∼ 15 kpc mainly constitutes the bulge component of the progenitor galaxy. Of course, the total bulge mass and progenitor orbit are uncertain, but such a bulge remnant should be detected by integral field spectroscopic and/or photometric observations around the predicted position. In addition, if high-velocity dispersion (probably induced by the MBH) occurred in the central region, the bulge component would be easily recognized (see e.g. Seth et al. 2014).

Figure 20.

Phase-space mass-density distributions (or position-velocity diagrams) of the progenitor's bulge and M31 stars. The inclination of the disc model is THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |⁠.

5.4 Extended stellar shell

Our unprecedented highly resolved simulation of the minor merger enables to predict a faint but huge stellar structure outside the western shell (see Fig. 17). The outer western shell (hereafter OWS) originates from outermost region of the initial progenitor galaxy. As shown in Fig. 17(b) and 19(c), the metallicity is lower in the OWS than in the GSS and both shells. On the other hand, the westernmost side of the GSS is sourced from the outermost region of the initial progenitor galaxy and appears as a broad GSS structure. The OWS metallicity will also limit the initial metallicity distribution of the satellite progenitor. Estimating the surface brightness of the OWS is important for further observations. A faint extent of metal-poor component has appeared in PAndAS observations (see fig. 2 of Martin et al. 2013), which may correspond to our simulated faint shell. This correspondence requires validation by additional spectroscopic observations of the faint component.

Fig. 21 shows the phase-space mass-density distributions in the western shell and OWS. The analysed area is | $180^\circ <\theta _a<230^\circ$ | with radii R > 0| $_{.}^{\circ}$ |5. This figure is constructed from the data set of the highest resolution run in the THICK model | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |⁠. To compare the phase-space distributions of the disc and Plummer progenitors, we use data with very similar mass resolutions in the two cases. Fig. 21(a) shows the western shell at | $R < 2^\circ$ | and the OWS at | $R < 3^\circ$ |⁠. Both shells are clearly distinguished by their phase-space mass-density distributions (Fig. 21b). Panels (c) and (d) of Fig. 21 reveal a similar structure in a spherical progenitor merger. The disc and Plummer models differ in their OWS phase-space distributions; specifically, the latter model exhibits a symmetric pattern in V_los (Fig. 21c) whereas the former shows an asymmetric pattern (Fig. 21a).

Figure 21.

Phase-space mass-density distributions of the simulated western shell and OWS assuming the disc (a and b) and Plummer (c and d) progenitors. Panels (a) and (c): line-of-sight velocity distributions. Panels (b) and (d): phase-space distribution centred on M31's centre. The inclination of the disc model is THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |⁠.

Fig. 22 shows the line-of-sight velocity and spatial distribution of the velocity dispersion in the observed frame. The line-of-sight velocity dispersion of σ_los ≃ 0 reveals clear edges of the western shell and OWS (Fig. 22a). The OWS particles have experienced two pericentric passages, as particles in the eastern shell. Therefore, on the phase-space distribution, the particles in the OWS and eastern shell exhibit the same phase.

Figure 22.

Distribution of line-of-sight velocity dispersion (a) and velocity (b) in the simulated western shell and OWS. Symbols and lines are those of Fig. 5(a). The inclination of the disc model is THICK7 with | $(\phi ,\theta )=(-15^\circ , 30^\circ )$ |⁠.

At the OWS, the disc and spherical models differ primarily by their distances from us. In the Plummer model, the OWS corresponds to a semicircular arc at the bottom-right of Fig. 8(f) (⁠| $-3^\circ <\xi <0^\circ$ | and −45 kpc < d_M31 < 30 kpc). In the disc model, most of the stellar components in the OWS spread out only in the foreground of M31 (⁠| $-3^\circ <\xi <0^\circ$ | and −40 kpc < d_M31 < 0 kpc).

We here summarize the extended stellar shell. The surface brightness of the extended shell is almost flat and requires a V-band detection limit 3–4 mag deeper than the apparent magnitude at the western shell on the minor axis of M31's disc. The shell is observed in both disc and spherical progenitor models. The morphologies of the progenitor galaxies in the two models differ by their distance between the OWS and us. The stars in the extended stellar shell have relatively low metallicities with small azimuthal variation, although the metallicity in the azimuthal direction varies widely in the simulated western shell. Future observations of the extended stellar shell (beyond the currently observed region) would further constrain the progenitor model.

5.5 Gas distribution

The observed structures in the GSS favour a rotating disc galaxy as the progenitor. As disc galaxies frequently contain gaseous components, we can predict that such components are stripped and dispersed in M31's halo. H i observations around M31's disc have revealed high-velocity H i clumps that aligned the GSS with an offset of ∼15 kpc and a similar line-of-sight velocity to the GSS (Westmeier, Braun & Thilker 2005; Westmeier, Brüns & Kerp 2008; Lewis et al. 2013). Only recently, another gaseous component in the GSS was detected in the absorption spectra of a background active galactic nuclei source (Koch et al. 2015). In addition, the gaseous rings of M31's disc could have been formed by a recent gaseous interaction (Gordon et al. 2006), reminiscent of the past gaseous interaction of a gas-rich progenitor. The origin may be revealed by hydrodynamical simulations.

6 SUMMARY

Through detailed simulations of the merger event and comparisons with observed data, we have strictly constrained the physical quantities of M31 and the infalling progenitor, including the gravitational potential of M31, the progenitor orbit and progenitor mass. However, the morphology (and dynamics) of the GSS progenitor galaxy has not been detailed here. By simply analysing the stellar count maps of the GSS in M31's halo, we characterized the asymmetric surface brightness profile across the GSS, which constrains the morphology of the GSS progenitor. We also perform the first large systematic survey of a minor merger with a disc satellite progenitor galaxy.

We identified the parameter space that properly reproduces the asymmetric surface brightness of the GSS on the plane (ϕ, θ), which defines the inclination angle of the initial disc. The structure was best reconstructed by the thick disc model (R_d = 1.1 kpc, Z_d = 0.52 kpc). The dynamically hot disc model cannot easily reproduce the eastern edge of the GSS, because the GSS is broadened by velocity dispersion in this model. On the other hand, the thin disc model struggles to reproduce the broad western structure, because it generates an excessive dynamically cold component.

Finally, we summarize our four predictions gained from the highly resolved simulations, which could be verified in future observations. First, the progenitor's bulge currently occupies the eastern shell and foreground of the disc of M31 and the two structures are distinguishable in the phase-space mass-density distributions. Secondly, we expect clear metallicity differences in the merger remnants, because the metallicity clearly differed in the azimuthal direction at approximately 3| $_{.}^{\circ}$ |5 < R < 4| $_{.}^{\circ}$ |5. The western shell also exhibited clear metallicity differences in the azimuthal direction. Thirdly, an extended stellar shell should reside outside the western shell. This extended shell should be detectable in photometric observations if the detection limit in the V band is 3–4 mag deeper than the apparent magnitude of the western shell on the minor axis of M31's disc. Finally, the western and extended shells contain clearly different stellar populations and observations of their metallicities and/or distances would further constrain the progenitor model.

Acknowledgments

We thank M. J. Irwin for allowing us to use the observed data. We would like to thank the anonymous referee greatly for fruitful and helpful suggestions to improve the paper. We would also like to thank the Center for Computational Sciences, University of Tsukuba, where all of the numerical calculations are carried out on. This work was supported by Grant–in–Aid for JSPS Fellows (T.K. 26.348), and the JSPS Grant–in–Aid for Scientific Research (A)(21244013) and (C)(25400222). This work was also supported by the Japan Science and Technology Agency's (JST) CREST programme entitled ‘Research and Development of Unified Environment on Accelerated Computing and Interconnection for Post-Petascale Era’. T. Kawaguchi was supported in part by a University Research Support Grant from the NAOJ.