Constraining Local Group Dark Matter Using M33’s Past Orbit A Thesis Presented by

Marion I. P. Dierickx to The Department of Astronomy in partial fulfillment of the requirements for the degree of Bachelor of Arts

April 2012 Harvard University

Constraining Local Group Dark Matter Using M33’s Past Orbit Marion I. P. Dierickx [email protected]

ABSTRACT

We compare results derived from a newly developed semi-analytic model of Local Group orbits with output from GADGET-2 simulations. Existing constraints on the quantity of intergalactic dark matter in the Local Group are limited to the classic timing argument approach. Here we use M33’s past orbit as a probe for the diffuse dark matter distribution in the Local Group. Given the parameter space granted by uncertainties in M31’s proper motion and diffuse intergalactic dark matter, our semianalytic model finds plausible orbits for the Milky Way, M31 and M33 galaxies over the past 5 Gyr. The system is treated as a three-body problem of particles moving on a background potential of superimposed bulge, disk and Navarro, Frenk, & White (1997) dark matter halo profiles. Using collisional and tidal stripping constraints, an upper limit of ∼ 1 Milky Way mass is found for the amount of diffuse dark matter present in the Local Group. Previous studies have found evidence for a close M31-M33 past encounter in the tidal disturbance of M33’s disk. Using our models we evaluate the effect of such encounters on the simulated evolution of the M33 globular cluster population. Our results indicate that at most ∼ 3 × 1010 solar masses were stripped off of M33 due to interaction with M31 over the past 5 Gyr. Orbital solutions are then computed and confirmed with GADGET-2 and yield observable dynamic signatures in the globular cluster record. We show that some of M31’s outer halo clusters may have been accreted from M33 and propose their identification as a diagnostic of past Local Group orbital history.

–2– Contents

1 INTRODUCTION

3

2 A SEMI-ANALYTIC MODEL OF THE LOCAL GROUP

8

2.1

Model architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2

Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.3

Parameter space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4

Model Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.5

Globular clusters as tracers of past accretion . . . . . . . . . . . . . . . . . . . . . . 19

2.6

Creating a globular cluster distribution

. . . . . . . . . . . . . . . . . . . . . . . . . 24

3 THE LOCAL GROUP WITH GADGET-2 3.1

GADGET-2 background

3.2

GADGET methods

8

30

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4 RESULTS

36

4.1

Successful semi-analytic runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4.2

Comparison of semi-analytic and GADGET results

5 DISCUSSION

. . . . . . . . . . . . . . . . . . 38

50

–3– 1.

INTRODUCTION

As our “cosmic backyard,” the Local Group of galaxies is the most readily observed laboratory of galactic formation and evolution. The Milky Way, M31 and their roughly thirty satellites form a relatively isolated system with little effect from the Hubble flow, with the nearest galaxies similar in mass over 3 Mpc away (Binney & Tremaine 1987). Indeed, M31 (also known as Andromeda) is known to be approaching the Milky Way at a speed of 116 km/s (Brunthaler et al. 2005). Given their current separation of 780 kpc (Sherwin, Loeb & O’Leary 2008), the resulting eventual merger is predicted to occur in ∼5 Gyr (Cox & Loeb 2008). In addition to this galactic merger taking place in our Local Group laboratory, an important accretion event is unfolding simultaneously. The third most massive member of the Local Group, M33 (also known as the Triangulum galaxy), is believed to be in the process of accretion by its larger neighbor M31. The current configuration of the three main Local Group galaxies is shown in Figure 1. Large spirals like the Milky Way and M31 are thought to form through continual accretion and mergers of smaller galaxies such as M33 (McConnachie et al. 2009). Vestiges of accreted satellite dwarf galaxies can be resolved and provide a test for theories of hierarchical structure formation. At the same time, mapping the Milky Way and characterizing stellar populations within stellar streams allows astronomers to probe past satellite accretion events. Yet even within the Local Group, many fundamental galactic properties remain poorly characterized. Hierarchical structure formation, for instance, predicts substantially more dark matter sub-haloes than the number of observed dwarf satellites quoted above; this conundrum has come to be known as “the missing satellite problem.” Galactic cartography is a field where large uncertainties remain, most notably about the structure of our own Milky Way. The prominent bar at the center of our Galaxy was confirmed and characterized as recently 2005 with the GLIMPSE survey (Churchwell & GLIMPSE Team 2005). Due to difficulties in the kinematic mapping of HI clouds, uncertainties about the number of spiral arms present in the Milky Way remain (Benjamin 2008), although important progress has recently been made thanks to Very Large Baseline Array (VLBA) measurements of masers (Brunthaler et al. 2011). Similar high-precision interferometry

–4–

Fig. 1.— Relative configuration of the three main Local Group members.

of two water masers by Brunthaler et al. (2005) have yielded precise proper-motion estimates for M33. However, such measurements have not yet been implemented for our “sister galaxy,” due to a lack of known maser sources until the recent detections by Darling (2011). A number of indirect methods have however been employed to constrain its transverse velocity components. Noting that M33’s calculated past orbit depends on M31’s 3-dimensional velocity vector, Loeb et al. (2005) make use of Brunthaler et al. (2005)’s M33 proper motion measurements to simulate the M31-M33 system backwards in time. Given that the radial motion of Local Group galaxies are well known from Doppler measurements, by requiring M33’s survival to tidal disruption, they find constraints on the two-dimensional transverse components of M31. In addition, four mutually consistent statistical analyses by Van der Marel & Guhathakurta (2008) yield M31 proper motion estimates with uncertainties on the order of 10 µas/yr. These uncertainties introduce significant ambiguities in our understanding of past galactic orbits, and hence in our knowledge of the Local Group’s history. In addition to M31’s velocity vector, in my junior thesis work M33’s past orbit was found to

–5– be sensitive to the Milky Way’s dark matter halo mass. Yet the distribution of dark matter in the Local Group remains poorly characterized. The shape of the Milky Way’s dark matter halo is an ongoing controversy, with different studies characterizing it as spherical, oblate or triaxial (eg Watkins, Evans & An (2010)). As recently as 2008, Xue et al. (2008) derived a very low Milky Way mass estimate of 1 × 1012 M based on SEGUE stellar kinematic samples. In contrast, Reid et al. (2009) used VLBA proper motion measurements to produce a radically higher estimate, with the Milky Way’s dark matter halo found to be similar in mass to M31’s (about 1.6 × 1012 M ). Evaluating dark matter halo mass relies on methods falling into two broad categories (Watkins, Evans & An 2010): intrinsic criteria such as finding the surface brightness profile of the stellar halo or the amplitude of the gas rotation curve; and methods relying on external tracers such as the kinematics of globular clusters and satellite galaxies, or the tidal radii of dwarf spheroidals. Thanks to gas rotation curves, both the Milky Way’s and its sister M31’s dark matter masses are reasonably well known within a few tens of kiloparsecs of their center (e.g. Rohlfs & Kreitschmann (1988)). By studying the kinematics of the satellite population in the outer halo, Watkins, Evans & An (2010) produce mass estimates for the Milky Way and M31 within 300 kpc from their center. Little is known however about Local Group intergalactic dark matter. So far, the only mass estimate for the entire Local Group comes from a classic approach called the timing argument, first presented by Kahn & Woltjer (1959). The authors note that M31 is approaching the Milky Way at ∼ 120 km/s and that these two galaxies’ masses dominate the Local Group baryonic matter at the ∼ 90% level, warranting treatment as an isolated system of two point masses on a radial orbit. Assuming that the Local Group constitutes a “physical unit,” the Milky Way and M31 must have formed closer together than they are now, and must hence be on an orbit around their center of gravity. The age of the Universe provides a natural upper bound for the period of this orbit, which through Kepler’s laws is translated into a lower bound of 1.8 × 1012 M for the total mass of the Local Group. This early estimate is flawed due to the use of M31 and Milky Way masses that were off by an order of magnitude, and a separation off by 25%. Unaware of the existence of dark matter, and noting that an implausibly high number of small stellar clusters would be needed to account for the large quantity of missing mass, Kahn & Woltjer (1959) attribute this intergalactic mass

–6– to gas. Today the intergalactic mass is believed to consist of dark matter. For a better constraint, we can revisit the timing argument following Binney & Tremaine (1987) with more accurate input observables. The Milky Way and M31 are roughly spherically symmetric mass distributions with a separation that is larger than their size, and they can hence be modeled as two isolated point masses. The following general parameterization of Kepler’s equations can be used to solve this two-body problem (Binney & Tremaine 1987):

r = a(1 − e cos η) r t=

a3 (η − e sin η) + C GM

(1)

(2)

Here r is the separation between the two galaxies as a function of time t, a is the semi-major axis of the orbit, e the eccentricity, η the eccentric anomaly, G the gravitational constant, M the summed mass of the two bodies and C a constant. We set C equal to zero (Binney & Tremaine 1987) by requiring r = 0 at t = 0, which means e = cos η(0) = 1, i.e. the orbit it radial. We need a third equation in order to solve for M (and a and η along the way): dr =v= dt

r

GM sin η a 1 − cos η

(3)

At the present time we know the separation r and the radial velocity v. An upper bound on the time t since the start of orbital motion is the age of the Universe, 13.7 Gyr, yielding a lower bound on the mass. We can solve this system of simultaneous equations numerically to obtain M ≥ 4.6 × 1012 M . Given the masses of the main galactic components given in Table 1, the orbital timing argument leads to an excess mass of ∼ 1.5 × 1012 M . According to this fairly crude approach, the M31 radial velocity and distance observed today imply a true Local Group mass higher by about a Milky Way mass than the summed masses of our Galaxy and M31. In addition to intergalactic dark matter, satellite galaxies may contribute to this predicted excess. M33, by far the largest Local Group satellite, has a mass of 2.2 × 1011 M (Seigar 2011), while most dwarf satellites, for example the Magellanic Clouds, have masses on the order of ∼ 109 , 1010 M (Kunkel,

–7– Demers & Irwin 1996). As a result, it seems that satellite galaxies in the Local Group can account for a fraction of the predicted additional mass, but not for the full excess. The aim of this work is to constrain the mass distribution of the diffuse intergalactic dark matter in the Local Group. Extending Watkins, Evans & An (2010)’s method of studying the kinematics of outer halo satellites, we use M33’s orbit as a visible probe of the underlying dark matter distribution. By examining their effect on M33’s past orbit, along the way we explore some of the consequences of unknowns in the intergalactic dark matter distribution and M31’s proper motion on Local Group hierarchy and dynamics. The stripping of globular clusters from M33 due to potential close passages near M31 can provide a visible tracer of M33’s past orbit. The methods used to predict globular cluster trails are based on preliminary work done for my junior thesis entitled “Modeling Globular Cluster Orbits under the Combined Gravitational Influence of Local Group Members.” The approach we take is two-pronged. First the Local Group is modeled semianalytically, allowing for computational flexibility and many runs in searching the parameter space. Successful parameter sets are then confirmed by simulating the orbits with the N-body GADGET code.

–8– 2.

A SEMI-ANALYTIC MODEL OF THE LOCAL GROUP

2.1.

Model architecture

Components The model is composed of four building blocks: the Milky Way, M31, Local Group dark matter, and M33. In 1997, Navarro, Frenk & White found that equilibrium dark matter halos could be fit with a universal density profile whose shape shows no sensitivity to halo mass. The NFW density profile is fully specified by its characteristic density δc , and has power-law dependencies on the radial distance r (Navarro, Frenk, & White 1997):

ρ(r) =

3H02 Ωm δc (1 + z)3 z 8πG Ωm cx(1 + cx)2

(4)

where z is the halo formation redshift, x = r/rvir , c is the halo concentration parameter, Ωm = 0.27 is the matter density, and the Hubble constant H0 is 71 km/s/Mpc (H0 = 100h, where h = 0.71) (following Loeb et al. (2005)). The concentration parameter is related to δc through the following expression: δc =

∆c c3 3 ln(1 + c) − c/(1 + c)

(5)

The overdensity ∆c is given in Barkana & Loeb (2001): ∆c = 18π 2 + 82d − 39d2

(6)

Ωm (1 + z)3 Ωm (1 + z)3 + ΩΛ + Ωk (1 + z)2

(7)

Here d = Ωzm − 1 and Ωzm =

Adopting standard Λ cold dark matter cosmological parameters, we take Ωk = −k/H02 = 0 in a flat Universe, and the cosmological constant density parameter ΩΛ = 1 − Ωm = 0.73. We take the

–9–

Fig. 2.— The scaled density of M33’s NFW profile (assuming a virial mass of 2.2 × 1011 M ) versus the logarithm of the radial distance scaled by the virial radius. This illustrates the double powerlaw profile of the NFW model. The cutoff radial distance between the two power laws is the scale radius, and here it is approximately indicated by the dashed line such that log(rscale /rvir ) ≈ −0.6. The concentration parameter is given by rvir /rscale , which gives cT ≈ 4, and so as expected we recover the value given in Table 1.

virialization redshifts for the three galaxies to be z = 1 (Loeb et al. 2005), while the Local Group collapse halo redshift is taken to be z = 0.5. This redshift value allows for integrations back in time up to 5 Gyr before the present. The NFW profile for M33 is shown in Figure 2, using the values for M33’s parameters listed in Table 1.

Having set up a framework for the formation of dark matter haloes, we can now proceed to determine the gravitational acceleration for an NFW profile. The characteristic parameter of a halo is its mass M , which yields the following virial radius for a halo collapsing at redshift z (Barkana & Loeb 2001):  rvir = 0.784

M 8 10 h−1 M

1/3 

Ωm ∆c Ωzm 18π 2

−1/3 

1+z 10

−1

h−1

kpc

(8)

– 10 – We use the circular velocity profile (Vc (r)) from Navarro, Frenk, & White (1997) to find the radial gravitational attraction for the NFW model:

Vc (r)2 = Vc2

1 ln(1 + cx) − cx/(1 + cx) x ln(1 + c) − c/(1 + c)

which gives gr = −Vc2 (r)/r

(9) (10)

where Vc is the halo’s circular velocity, given by Vc2 = GM/rvir (Barkana & Loeb 2001). So finally, we combine these to obtain the gravitational acceleration from an NFW dark matter halo:

gr = −

GM ln(1 + cx) − cx/(1 + cx) r2 ln(1 + c) − c/(1 + c)

(11)

In addition, the disk and bulge components of M31’s potential are modeled using the following forms (Xue et al. 2008; Shattow & Loeb 2009):

Exponential disk: Hernquist bulge:

GMA,disk (1 − e−r/bA ) r GMA,bulge Φbulge (r) = − r + c0,A Φdisk (r) = −

(12) (13)

The accelerations are then simply given by the negative gradient of the potential. Bulge/disk components for the Milky Way are omitted, as M33’s orbit is most affected by M31 and any disk and/or bulge contribution for the Milky Way is expected to be very small. The values used for the different model components are summarized in Table 1.

– 11 –

Table 1. Parameter

List of parameters used and their values.

Meaning

Value

ΛCDM cosmological parameters (Loeb et al. 2005) H0 Hubble constant 71 km/s/Mpc Ωm Dark matter density 0.27 M31 parameters (Sherwin, Loeb & O’Leary (2008), Shattow & Loeb (2009)) MA Virial mass 1.6 × 1012 M MA,bulge Bulge mass 2 × 1010 M MA,disk Disk mass 8 × 1010 M bA,bulge Bulge scalelength 1.8 kpc bA,disk Disk scalelength 5.5 kpc c0,A Disk scaleheight 1.14 kpc c Halo concentration 12 Milky Way parameters (Shattow & Loeb 2009) MM Virial mass c Halo concentration

1.5 × 1012 M 12

M33 parameters (Seigar (2011), Ciardullo et al. (2004)) MT Virial mass 2.2 × 1011 M MT,bulge Bulge mass 1.14 × 108 M MT,disk Disk mass 3.81 × 109 M bT,bulge Bulge effective radius 0.39 kpc bT,disk Disk scalelength 1.7 kpc c0,T Disk scaleheight 0.175 kpc

– 12 – Overall procedure Each run is initialized at t = 0 with M31, the Local Group and the Milky Way having density profiles as explained above. The Local Group halo is centered at the current Local Group barycenter (neglecting any satellite galaxies less massive than M33 in this calculation), as the highest concentrations of dark matter would naturally be expected at the system’s center of mass. To facilitate computation, the Local Group contributes a static background potential anchored at the center of mass throughout the integration. This is a valid approximation because the halo is massive and starts out at the center of mass with zero velocity, so its later motion is expected to be small. This was verified with GADGET in Section 4.2. Given these six-dimensional phase space initial conditions, we integrate back in time to t = −5 Gyr. The galaxies are treated as point particles moving in the potential fields produced by summing NFW profiles as well as the bulge and disk components or M31. We have the following simple system of equations of motion: 







d  ~v   ~a   =  dt ~x ~v

(14)

This system is integrated with a fourth-order Runge-Kutta method. After checking the convergence of a number of different timesteps, we settle on using a timestep of 1 Myr as it is sufficiently accurate and allows for efficient computing. Table 1—Continued Parameter

Meaning

Value

cT

Halo concentration

4

Note. — In this paper the subscripts A, M and T refer to M31/Andromeda-, Milky Way- and M33/Triangulum-related quantities respectively.

– 13 –

t=0

Track: - M33 tidal stripping - GC coordinates

Orbital integration

t = - 5 Gyr

Model predictions: - M33 mass and radius - GC trail

Current observables & parameters

Record phase-space coordinates

Initial conditions

Assign M33 mass and GC distribution

Fig. 3.— General procedure for a run of the analytical model.

At t = −5 Gyr, the phase-space coordinates of all model components are recorded. While its orbit is still computed as a point particle, M33 is assigned a virial mass and a corresponding NFW profile. This allows tracking the tidal stripping of M33 as it passes near M31. Since M33’s orbit is that of a point particle (with no dynamical friction) the trajectories are reversible and at the end of the run, back at t = 0, we obtain the initial phase-space conditions. We then evaluate the extent of M33’s stripping by checking whether the end-mass lies within the error bars of our value for MT , (2.2 ± 0.1) × 1011 M . As an additional constraint on successful models, we use a no-collision condition (see Section 2.4). If the run is successful at predicting appropriate tidal stripping for M33, we go back to t = −5 Gyr and assign a virialized population of 200 globular clusters following an NFW profile. The methods for constucting this population are outlined in Section 2.2. The overall structure of a successful run is summarized in Figure 3.

– 14 – 2.2.

Initial conditions

Throughout the integration we use a Cartesian coordinate system centered at the initial Local Group barycenter. The Local Group dark matter halo is hence held fixed at (0, 0, 0). Current observables are taken for the phase-space coordinates used to initialize each run (see Table 2). Distances, galactic longitude and latitude are converted to galactic Cartesian coordinates by adapting standard spherical to Cartesian conversion, as shown in equation 15.

x = d cos b cos l,

y = d cos b sin l,

z = d sin b

(15)

The location of the Local Group’s barycenter is then computed, and the origin translated to that location. Matrix conversions from proper motions and radial velocity to galactic space-velocity components [U, V, W ] are available in Johnson & Soderblom (1987). The input for M31 is velocity in the north and west direction, which we convert to proper motion using the standard δ˙ = µN = vN /D and α˙ = −µW / cos(δ) = −vW /(D cos δ) relations. Once the velocity components are known in the galactic Cartesian frame, we correct for the Local Standard of Rest (LSR) motion, which following Brunthaler et al. (2005) is taken to be 236 km/s. A correction for the solar peculiar motion of [10, 5.25, 7.17] km/s (Brunthaler et al. 2005) is also applied. Next we can proceed to calculate the velocity vector of the Local Group barycenter. We switch to the Local Group center of mass (CM) frame by vector subtraction, ~r = ~rgalactic − ~rCM , and so the Local Group dark matter profile is at rest in this frame.

2.3.

Parameter space

We defined a grid in the parameter space as input to the analytic model. The four free parameters are: • The Local Group mass. In practice it is defined as a mass parameter µLG multiplying the summed observed masses of the three main Local Group galaxies: the Milky Way, M31, and

– 15 – M33. The range of µLG in the first set of runs is from 0 (no Local Group diffuse dark matter) to 2. This upper limit is generous, as values of µLG larger than approximately 1.5 yield trial galactic trajectories that are very disordered and include collisions. • M33’s original mass. In this work the standard value for M33’s virial mass is taken to be MT = (2.2 ± 0.1) × 1011 from Seigar (2011). However, M33’s starting mass, the mass which is attributed to the galaxy after 5 Gyr of integration back in time, is a free parameter. Here again we use a mass parameter that multiplies the standard value of MT : values higher than approximately 2 do not produce sufficient stripping for M31’s final mass to be within the error bars at the end of the run. • M31’s proper motion in the western and northern directions, vW and vN . M31’s heliocentric radial velocity is well-known due to Doppler measurements, and is taken to be −301 km/s (Courteau & van den Bergh 1999). However, the fact that its transverse motion components remain poorly constrained causes major uncertainties in the dynamical past and future of the Local Group (Loeb et al. 2005). Precise measurements of M33’s proper motion have been made with VLBI (Brunthaler et al. 2005). These were used by Loeb et al. (2005) to place limits on M31’s proper motion based on the constraint that M33’s thin stellar disk be minimally disrupted. Van der Marel & Guhathakurta (2008) present four independent, mutually consistent statistical methods to determine M31’s velocity vector. They yield an average of hvW i = −78 ± 41 km/s and hvN i = −38 ± 34 km/s. We base our parameter space for M31’s proper motion on these error bars, and sample the allowed velocity range at fixed intervals. Darling (2011) recently reported the discovery of five water masers in M31 with the Green Bank telescope, which will soon allow for precise proper motion measurements. The motion of M31*, the compact bright radio source at M31’s nucleus, can now be measured with VLBI and will provide more direct proper motion estimates than masers (Mark Reid, personal communication).

Under the assumption that each parameter’s effect on the orbits varies smoothly, an initial

– 16 –

M31 M33

R

rt v

Fig. 4.— Illustration for the calculation of the tidal radius rt .

4 × 5 grid was defined to sample the parameter space, as shown in Table 3. The successful runs found though this grid each are representative of a small class of models, which then warrants more detailed sampling. This grid represents an initial set of 54 runs, yielding a computation time of a little over two hours. In total two sets of runs were conducted, the second set with refined grid spacing, as can be seen in Table 4.

2.4.

Model Constraints

Tidal stripping As M33 moves on its trajectory near M31, it experiences tidal forces across its halo. To derive the condition for tidal disruption, we consider the change in force across M33 as compared to its binding force. A parcel of mass in M33 situated a distance r from its center will experience the

– 17 –

Table 2.

Observables used as initial coordinates.

Parameter

M31

M33

l (degrees) b (degrees) Distance (kpc) Radial velocity (km/s) α˙ (µas/yr) δ˙ (µas/yr)

121.174 (NED) -21.573 (NED) 785 ± 25 (McConnachie et al. 2005) -301 (Courteau & van den Bergh 1999) Free Free

133.61 (NED) -31.331 (NED) 794 ± 23 (McConnachie et al. 2004) -179 (Corbelli & Schneider 1997) 23.2 (Brunthaler et al. 2005) 7.5 (Brunthaler et al. 2005)

Note. — In row 3, the distance estimates are based on Tip of the Red Giant Branch (TGRB) analysis, in good agreement with values from Cepheid variables (McConnachie et al. 2004). In rows 5 and 6, the proper motion components of M31 are left as free parameters.

Table 3.

List of parameters and their values.

µLG

µT

vW (km/s)

vN (km/s)

0.0 0.5 1.0 1.5 2.0

1.00 1.25 1.50 1.75 2.00

-119 -98.5 -78.0 -57.5 -37.0

-73.0 -56.0 -38.0 -21.0 -4.0

Note. — Parameter grid for initial set of 54 = 625 runs. µLG is the mass parameter multiplying the summed observed masses of the three main Local Group galaxies to give a Local Group halo mass. µT is the mass parameter multiplying the standard value of MT to give an initial M33 mass. vW and vN are the present-day M31 proper motions in the western and northern directions. With each run taking roughly 12 seconds, this grid represents a total computing time of a little over two hours.

– 18 – following gravitational forces:  GMA (< R) r R2 R   GMA (< R) r R2 R  Bound: Unbound:

< >

GMT (< r) r2 GMT (< r) r2

(16) (17)

R is the separation between the two galaxies’ centers, MA (< R) the M31 mass interior to radius R, MT (< r) the M33 mass interior to the orbital radius r of the mass parcel under consideration. These terms are illustrated in Figure 4. The tidal radius is then given by equalling the two forces:  rt =

MT (< rt ) MA (< R)

1/3 R

To find the mass inside of a given radius, we use equation 10, where gr =

(18) GM (
and so in general

M (< R) = r · Vc2 (r)/G. Due to the dependence on rt embedded in the MT (< rt ) term, equation 18 can only be solved numerically. This is done using a modified Newton-Raphson routine in IDL. At each timestep we compute the tidal radius and update it to the lowest value calculated so far, as material that is stripped is lost, and as a result the tidal radius can never increase. The final value of the tidal radius is recorded, and the mass inside is computed. We obtain a constraint by requiring that this mass lie within the error bars of M33’s mass estimate from Seigar (2011).

No-collision conditions Sets of parameters leading to past collisions between any of the three Local Group galaxies are eliminated, as there is no observational evidence for such events in the past 5 Gyr. A collision between galaxies i and j is defined as the virial haloes of both galaxies involved overlapping substantially: ri,j ≤

rvir,i + rvir,j 2

(19)

In the third set of runs, an analogous no-collision condition was implemented for the Local Group halo and M33. Since the tidal stripping from Local Group dark matter is not computed, in previous

– 19 – runs trajectories where M33 has a close passage near the Local Group halo, leading to M33’s destruction, were not being eliminated. Further development of this model will involve including tidal stripping of M33 due to the Local Group halo.

2.5.

Globular clusters as tracers of past accretion

Lessons from the Milky Way globular cluster population Giant spiral galaxies like the Milky Way or M31 are believed to harbor a few hundred globular clusters (GCs), which are concentrated spherical groups of bound Population II stars. In their catalog of M33 globular clusters, Sarajedini & Mancone (2007) identify 255 confirmed by the Hubble Space Telescope (HST) and high-resolution ground-based observations. By making use of HST’s high angular resolution, almost 500 globulars have been cataloged for M31 (Barmby & Huchra 2001). Fewer are known for the Milky Way; as of June 2011, 157 GCs have been found in our galaxy (Harris online catalog, Harris (1996)). It is speculated that 10-20 more have yet to be found due to obscuration from the galactic bulge. As GC stars are believed to have formed simultaneously, they share chemical abundances and enrichment histories, and have hence been valuable in the fields of galactic archaeology and stellar evolution (Harris 1996). Following the work done in the context of my junior thesis in Astronomy 98, I hope to use tidally stripped or accreted globular clusters as observable tracers of M33’s past orbit. Using globular cluster trails as a probe of past dynamics, we must be able to differentiate between globulars that formed in M31 in situ and globulars that have been accreted. Here we investigate the feasibility of this approach based on studies of what is known about the Milky Way globular cluster population. Recent comprehensive studies have attempted to estimate the number of accreted dwarf satellite galaxies based on the number of accreted GCs (Forbes & Bridges 2010; Mackey & Gilmore 2004). Mackey & Gilmore (2004) investigate sub-samples of the Milky Way GC population corresponding to the so-called “old halo,” “young halo,” and “bulk/disk” subsystems. Spatial and age distributions, age-metallicity relationships and orbital parameters are used to compare these

– 20 – three populations with GCs located in satellite dwarf spheroidal galaxies such as the Magellanic Clouds, Formax and Sagittarius. It has been proposed that “young halo” clusters in fact formed in such external galaxies that were later accreted. Mackey & Gilmore (2004) support this hypothesis based on their observation that external GCs are virtually indistinguishable from young halo clusters in terms of their horizontal-branch (HB) morphology, derived from Hubble Space Telescope photometry. HB morphology describes the proportion of stars found at different spectral types in the HB turnoff on the GC’s Hertzprung-Russel diagram, and is related to chemical composition as well as age; young halo objects tend to have redder HB morphologies. As another piece of evidence that “young halo” GCs originate in satellite galaxies, Mackey & Gilmore (2004) find that the core radius distribution of external clusters closely matches that of young halo GCs. In addition, young halo GCs extend to very large galactocentric radii, are found on higly eccentric and retrograde orbits, and possess similar frequencies of RR Lyrae stars as external clusters. Mackey & Gilmore (2004) conclude that roughly 40 objects in their sample are of extra-Galactic origin, most of them members of the young halo population. Forbes & Bridges (2010) build on the methods proposed in Mackey & Gilmore (2004) to further analyze the Milky Way GC population as compared to neighboring dwarf spheroidals. They look at the clusters’ age-metallicity relationship (AMR) in order to use their chemical enrichment history as an additional distinguishing factor. They find two diverging tracks in their AMR analysis; one that has a near-constant old age of ≈ 12.8 Gyr over the full metallicity range, and a bifurcation track to younger ages associated with the same metallicities. From the analysis in Forbes & Bridges (2010), it appears that a GC having a young age at a given metallicity is an indicator of possible accretion origin. From HB morphology models, Yoon & Lee (2002) had previously found evidence that the properties of accreted globular clusters of metallicity similar to that of a given group of in situ GCs are well-explained if the accreted GCs are younger. The authors offer an explanation for the Oosterhoff dichotomy, a long-known division of Galactic GCs based on the mean period of RR Lyrae variable stars. This dichotomy is also traceable in systematic differences in metal abundances

– 21 – and kinematic properties. Yoon & Lee (2002) suggest that the Oosterhoff dichotomy and related discrepancies in observed GC properties are explained by the in situ formation of GCs in the proto-Galaxy, combined with the later accretion of GCs from satellite systems. Reviewing these different studies, there seem to be at least four GC properties that provide evidence of either in situ or external origins. Below I summarize some of the observations that allow one to distinguish between these different GC sub-populations: • Age-metallicity relation. The location of globular clusters on the [Fe/H] vs. age graph reflects different chemical enrichment histories. Following Forbes & Bridges (2010), it appears that in situ GCs tend to be older and that lower metallicities may indicate accretion from low-mass galaxies. Thus, a GC’s young age for a given metallicity may indicate an accretion origin. • Horizontal branch morphologies. An HB morphology similar to those found in external GCs indicates a plausible accretion origin (Mackey & Gilmore 2004; Forbes & Bridges 2010). • Core radii distribution. Mackey & Gilmore (2004) find very distinct core-radius distributions for the young halo, old halo and bulge/disk GC subpopulations in the Milky Way. Core radii in the young halo group are found to be significantly more extended, matching the observed distribution for external GCs. In addition, very extended GCs are found in the Milky Way and not in the external sample of dwarf satellites, raising the possibility that these clusters are being disrupted following disturbance during accretion. • Orbital kinematics. The nature of globular cluster orbits may be the most definite diagnostic parameter for a captured origin (Yoon & Lee 2002; Mackey & Gilmore 2004; Forbes & Bridges 2010). Specifically, retrograde, high-eccentricity and/or high-energy orbits are good indicators of provenance from a satellite galaxy on a retrograde orbit before and during accretion.

– 22 – M31’s globular clusters Can the observations necessary to characterize these properties can be carried out for globular clusters in M31? The fact that M31 is roughly 8 times more distant than the most distant Galactic GCs is an obvious challenge and an important one, as two of the methods outlined above rely on observing individual stars within the cluster. To compute the AMR, the age is determined from looking at the main-sequence turnoff of the stars in the cluster. The color-magnitude diagram can also give the metallicity (Rich et al. 2005). At present, the necessary resolving of individual stars in the cluster can best be done with HST or with very large ground-based telescopes equipped with effective adaptive optics systems (Rich et al. 2005). In spite of these observational challenges, M31’s globular cluster population has been characterized in some detail. A spectroscopic study of several hundreds of M31’s GCs, conducted with the William Herschel Telescope, adds to the mounting evidence in favor of a bimodal distribution in metallicity (Perrett et al. 2002). In addition, these two GC populations are found to differ in spatial distribution and kinematics. It is possible that the metal-poor, less spatially concentrated population originates from accretion. In-detail photometry with HST allows Rich et al. (2005) to construct high-quality colormagnitude diagrams for a sample of M31 GCs. These clusters are found to be quite similar to the Milky Way population. A number of GCs are found with metallicity-HB type relationships resembling the HB morphology of Formax globular clusters. Similarly to Yoon & Lee (2002), the authors propose that this offset could be accounted for by an age difference, but they do not go as far as to suggest that this age offset may be due to the incorporation of younger GCs from accreting dwarfs. However, their findings are limited by the small sample size of their study, which includes only 19 objects. Mackey et al. (2010) take a different approach and used newly-released data from the PanM31 Archaeological Survey (PAndAS) to map a new sample of outer halo globular clusters in M31. Using Monte-Carlo methods, they analyze the likelihood of chance association of observed stellar stream structure with the GCs’ positions. Their study provides evidence that GCs in M31 are anisotropically distributed and preferentially aligned with underlying tidal debris features. A large

– 23 – majority of the observed GCs coincide with prominent substructure. These important results imply that a significant portion of the outer halo GCs have been accreted from satellite parent galaxies. In addition, spectra and/or color-magnitude diagrams were available for a small number of GCs in the new PAndAS sample. As a preliminary result, most of these objects were shown to be metal-poor, further supporting a possible accretion origin.

M33’s globular clusters Exploring the galactic outskirts of M33 with the Isaac Newton Wide-Field Camera, Huxor et al. (2009) discovered four remote GCs, that is, objects with radial distances to the center of M33 larger than approximately 10 kpc. They note that the three outermost clusters are located on the far side of M33 with respect to M31, suggesting that tidal interactions could have affected M33’s globular cluster distribution at large radii. Given the small number of GCs found, this observation is limited by small number statistics. PAndAS was used to search for new clusters in M33 as well as in M31. A wide-sky coverage allowed Cockroft et al. (2011) to look for outer halo GCs, reaching to a projected radius of 50 kpc. Interestingly, only one new GC was found in addition to the five previously known. The authors suggest that M33’s outer halo clusters were stripped off in a previous dynamical encounter with M31. To summarize, M33’s globular cluster population has been under-studied as compared to the Milky Way’s or the Magellanic Clouds’. While it is currently challenging to observe M31’s and M33’s GCs, present observing facilities have the capacity to carry out more detailed studies, and will likely do so in the near future to follow up on recent breakthroughs from PAndAS. First results indicate that 1) M31’s GC population is very likely to have a significant accreted component, the exact origin of this component (accreted dwarf spheroidals, or close encounter with M33) having yet to be determined; 2) it is possible that M33’s population has been influenced by tidal interaction with M31 in the past. Improved photometric mapping and spectroscopic studies are needed to disentangle these populations.

– 24 –

2.6.

Creating a globular cluster distribution

The methods outlined in this section are very similar to the work done in the context of my junior thesis, but have been updated in light of a few inaccuracies. Drawing from Section 2.5, at t = −5 Gyr we assign to M33 a population of 200 GCs following the NFW profile corresponding to the mass parameter µT (see Table 3). A realistic set of initial phase-space coordinates is needed for each object.

Position Spherical coordinates are sampled in order to match the NFW density profile in equation 4. Since the profile is spherically symmetric and there is a radial dependence only, φ and cos θ are sampled from a uniform distribution. We use acceptance-rejection sampling, a Monte-Carlo technique, to sample r from the probability density function associated with the NFW profile. This distribution is given by: P (r) =

1 4πr2 ρN F W (r) Mtot

(20)

where Mtot is the total mass, and the 4πr2 factor accounts for the P dependence on φ and θ. This function is plotted in the top panel of Figure 5. The rejection method functions as follows. First a sample r smaller than the virial radius of M33 is drawn from a uniform distribution. This maximal radius rv,T is generous. Indeed, it has been found that very few clusters in M33 are at radii larger than 50 kpc (see Section 2.5 as well as Cockroft et al. (2011)). For comparison, a virial mass of 2.2 × 1011 M corresponds to a virial radius of 94 kpc. For many runs, the virial radius will be larger due to larger values of the parameter µT . Next, an associated probability is drawn uniformly between 0 and the maximum of the probability density function. If this probability lies under the curve, then the associated radius is kept. This algorithm is repeated until enough points are found;

– 25 –

Fig. 5.— The top panel shows the probability density function as a function of radius for the globular clusters in M33, assuming MT = 2.2 × 1011 M . The bottom panel shows the radial distribution of a sample set of 200 GCs derived according to the probability density function above.

the resulting radial distribution is shown in the bottom panel of Figure 5. Velocity Here we discuss how to determine isotropic velocity vectors for each GC in our sample using the Eddington formula. There is no angular dependence, and so cos(θv ) and φv are sampled from a uniform distribution. Velocity magnitudes are to be sampled according to a probability

– 26 – density function derived from the globulars’ relative energy distribution function. Following Binney & Tremaine (1987), we can determine this function from each globular’s location within the spherically symmetric M33 NFW density profile. The relative energy and relative potential are defined as follows: 1 E = Ψ − v2 2

where

Ψ(r) = −ΦN F W (r)

(21)

For each radius r sampled above, these relationships allow for one-to-one mapping from a velocity v to a relative energy E. In 1916, Eddington derived the energy probability density function for a phase-space density as a function of radius only (Eddington 1916), as is the case here. In this discussion we will follow the derivation of Eddington’s equation found in Binney & Tremaine (1987). The density profile can be written as an integral of the probability density function f over the R velocity: ρ(r) = f d3~v . Switching variables to E, we obtain: Z ρ(r) = 4π

Ψ

f (E)

p 2(Ψ − E)dE

(22)

0

Since Ψ is a monotonic function of r, we can now rewrite ρ(r) as ρ(Ψ). This allows both sides of equation (18) to be differentiated with respect to Ψ, producing Abel’s integral equation on the right-hand side. The known solution to this equation is separated through integration by parts, and we obtain Eddington’s equation (Binney & Tremaine 1987): 1 f (E) = √ 8π 2

Z 0

E

1 d2 ρ dΨ √ +√ 2 dΨ E − Ψ E



dρ dΨ



 (23) Ψ=0

As we expect the NFW density to change very little at high r, where Ψ ≈ 0, the second term in this sum is zero. The ρ double derivative in the first term can be rewritten in a lengthier but more convenient way, such that it matches quantities that can be computed: d2 ρ dr d = 2 dΨ dΨ dr



dr dρ dΨ dr



dr d = dΨ dr

dρ 1 dr dΨ dr

! (24)

– 27 – Because the acceptance-rejection method does not require the probabilities to be normalized, all constants are ignored. After canceling the dΨ terms when combining equations 23 and 24, we obtain a practical version of Eddington’s formula: Z

r(Ψ=E)

f (E) = r(Ψ=0)

Working our way outwards, first

dΨ dr

and

dρ dr

1

d p E − Ψ(r) dr

dρ 1 dr dΨ dr

! dr

(25)

can easily be computed from equations 21 and 4. The

full derivative inside the integrand can therefore be done analytically. We now look at the bounds of the integral. Through equation 21, any given r sampled previously corresponds to a range of relative energies, determined by the range of allowed velocities at that radius, 0 ≤ v ≤ vesc (r). vesc (r) is p the escape velocity and is simply given by −2Φ(r). Looking for the r range that corresponds to this range in relative energy, we find that the equation that defines the upper bound of the integral, r(Ψ = E), is not invertible. The upper bound radius is hence solved for numerically using a modified Newton-Raphson method. The lower bound is where Ψ goes to zero, which occurs at infinite radii. The integral can now be computed for every value of E in the allowed range with an adaptive integration routine. This energy function then yields the probability density function for the velocity: f (v) = v 2 f (E(v)). This function and the resulting distribution are shown in Figure 6. As our set of globular clusters follows the dark matter distribution, we expect it to be selfgravitating and hence satisfy the virial theorem. This was verified by summing up the kinetic and potential energies of all globulars. The newly created spherical coordinates are in M33’s reference frame, with the center of the dark matter profile as the origin. The synthetic GC coordinates are then converted to our barycentric Local Group frame by vector addition.

Projection on the sky The final positions of the globular clusters are projected onto the sky as seen from Earth. This is done by inverting the coordinate transformations done in Section 2.2. We aim to produce an equal-area projection, where areas of the sphere of the sky are mapped onto a plane and scaled by

– 28 –

Fig. 6.— The top panel shows the velocity magnitude probability density function for the globular clusters in M33, again assuming MT = 2.2 × 1011 M . The bottom panel shows the velocity distribution of the same sample set of 200 GCs derived according to the probability density function above.

a constant scaling factor. This will ensure that we visually preserve the structure of possible tidal streams on the map. The Hammer-Aitoff projection, also known simply as Hammer projection, is based on Aitoff’s work on modified azimuthal maps and fits the entire sphere onto the plane by doubling longitudinal values from the central meridian. Unlike the Aitoff projection, the Hammer

– 29 – projection is also equal-area. (x, y) coordinates on the plane are given by:

x = y =

√ 2 2 cos(φ) sin(λ/2) p 1 + cos(φ) cos(λ/2) √ 2 sin(φ) p 1 + cos(φ) cos(λ/2)

where (λ, φ) are, respectively, the longitude from the central meridian and the latitude.

(26) (27)

– 30 – 3.

THE LOCAL GROUP WITH GADGET-2

While it has the advantage of offering flexibility and efficiency to rapidly obtain first-order approximations of the parameter space of interest, the model described above requires assumptions about the background galactic potentials. More importantly, our semi-analytic model is limited by a simplified treatment of galactic orbits and does not take into account dynamical friction or tidal effects between all bodies. These limitations can be alleviated by using an N-body simulation, which more faithfully renders gravitational and dynamical effects across each halo. In N-body simulations, Newton’s equations are solved and tracked for a large number of particles under their own self-gravity (Springel, Yoshida & White 2001). In the second phase of the project, we use GADGET-2 (Springel 2005), a freely available N-body code, to strengthen our parameterization of M33’s past interaction with the Local Group.

3.1.

GADGET-2 background

Increasingly powerful computing paired with the development of more efficient and sophisticated algorithms are fueling the field of numerical modeling. N-body simulations have become invaluable tools in theoretical cosmology and galactic evolution. Indeed, due to highly non-linear regimes and complex 3-dimensional geometries, direct numerical modeling is often the only approach to accurately representing gravitational and hydrodynamical effects (Springel 2005). Particle simulation methods have evolved greatly over the past few decades. In the 1970s, early N-body codes most often implemented direct summation methods, also known as particleparticle methods. In such schemes, forces from all particles are computed and accumulated for each particle. The equations of motion are then separated and solved through an integration scheme. In effect, the semi-analytic model described in Section 2 is a 4-body rendition of this method. The particle-particle method comes at a large computational cost for large N because the number of

– 31 – operations required to solve a system of N particles is of order N 2 . This approach is valuable however in specific cases such as collisional stellar dynamical systems (Springel, Yoshida & White 2001). Particle-Mesh methods (PM) were developed to improve computational speed. Here, forces are approximated on a mesh and treated as field quantities. Then the specific potentials and forces at each particle’s position are obtained by interpolating mesh-defined values. Solving the Poisson equation for the potential field can be done through Fourier transformations (Springel, Yoshida & White 2001). Here the number of required operations scales with N + Ng log Ng , where Ng is the number of grid points. PM methods have difficulty handling important density fluctuations, and as such their resolution is limited. One way this can be remedied is by using an adaptive mesh, where the grid size is refined at highly clustered regions. In the so-called Particle-Particle/Particle-Mesh (P3M) methods, the PM method is supplemented by separating long-range forces from short-range forces, which on scales below the mesh size are computed by direct summation (Springel, Yoshida & White 2001). Tree codes were pioneered in the 1980s as an alternative to PM methods (Springel, Yoshida & White 2001). Conceptually, they utilize the same force superposition scheme as P3M methods in that short-range and long-range forces are computed separately and superimposed. Tree algorithms divide the total particle volume into a hierarchical tree. Going down the tree when computing the force on a particle, nearby particles can be treated individually, while distant particles are grouped into massive pseudo-particles and their effect is approximated by computing their multipole moments (Springel, Yoshida & White 2001). Such tree algorithms allow for computational costs of order N log N and are intrinsically much more adaptive to high clustering than PM methods (Springel, Yoshida & White 2001). However, they are more costly in terms of storage requirements. Springel, Yoshida & White (2001) released the first serial and parallel versions of GADGET (GAlaxies with Dark matter and Gas intEracT) in 2001. Written in C, the software tracks the evolution of N self-gravitating, collisionless bodies and can optionally include gas dynamics. To compute gravitational forces, GADGET implements the standard Barnes and Hut tree construction,

– 32 – where the computational domain is divided recursively into three-dimensional cubes called oct-trees (because each parent cube contains eight sibling cubes) (Springel, Yoshida & White 2001). In regions of high density, particles may come close enough to each other to cause a divergence of their mutual gravitational force. As a result, the gravitational field is softened through a spline kernel method, which has the advantage that the force becomes exactly Newtonian for distances larger than a specified multiple of the softening length (2.8 in our case) (Springel, Yoshida & White 2001). GADGET allows different softening lengths for different types of simulated particles (dark matter vs. stellar matter for example). Fluids on the other hand are treated by means of Smoothed Particle Hydrodynamics (SPH). SPH is a particle-based, Lagrangian representation of fluid flow. Hydrodynamic forces and changes in internal energy are computed through assigning smoothing radii to discrete elements referred to as “fluid particles.” Properties of the fluid particle are smoothed over the smoothing length by kernel functions, which taper down to zero for distant elements. In this way the properties of each particle are influenced most heavily by its neighbors, resulting in a neighbor search similar in structure to the Barnes and Hut oct-tree walk (Springel, Yoshida & White 2001). The time integration for each particle is fully adaptive in that each particle carries its own, adaptive time step. In addition, GADGET featured a new parallelization strategy, implementing individual particle timesteps on distributed memory, massively parallel computers for the first time (Springel, Yoshida & White 2001). Unlike SPH fluid particles, gravitational particles interact with every other particle at every timestep, resulting in the need for efficient communication between processors. To this end GADGET utilizes the Message Passing Interface (MPI), a flexible communication scheme portable on a large variety of platforms. By far the most well-known computation carried out with GADGET is the Springel et al. (2005) Millenium Simulation, run in 2005 by the Virgo Consortium to test the cold dark matter model for hierarchical structure formation. This high-resolution simulation was a breakthrough in that it followed ten times as many particles as the previously largest N-body computation, and offered much improved time and spatial resolution (the Millenium Simulation tracked ∼ 1010 particles in a cube of side 500h−1 Mpc). Starting at redshift z = 127, the simulation and its post-processing

– 33 – identified key evolutionary links of structure formation and further consolidated standard ΛCDM cosmology (Springel et al. 2005).

3.2.

GADGET methods

In addition to the three galaxies whose orbits we integrate, we also require a dark matter halo representing the intergalactic Local Group dark matter. As a result, we seek to model four objects and existing programs to generate GADGET initial conditions for two objects had to be modified accordingly. Trial runs were conducted with dark matter haloes only; stellar matter in the form of bulge and disk components for M33 was input later for satisfactory orbits. Generating the galaxies Given a choice of resolution, the first step is to create the individual objects to be represented in the simulation. GADGET’s default dark matter halo profile is the Hernquist profile, with ρ(r) ∝ (r/rs )−1 (1 + r/rs )−3 , where rs is a scale radius (Hernquist 1990). The NFW profile used for the semi-analytic runs has the form ρ(r) ∝ (r/rs )−1 (1+r/rs )−2 . Unlike NFW, the Hernquist profile has the advantage that it does not diverge in total mass and thus does not have to be truncated at large radii. Given their similar form, NFW and Hernquist profiles have been shown to fit observed dark matter haloes equally well (e.g. McLaughlin (1999)). We therefore choose to stick with the GADGET built-in Hernquist profiles rather than attempting to modify the code structure. The profiles match each other closely in the inner regions of the halo (for small r); they do however deviate by a small amount in the outer regions (beyond ∼ 100 kpc). Here again the profiles are determined by their concentration parameter (see Table 1) and the circular velocity at the virial radius, V200 . As in the semi-analytic work, the Local Group dark matter is rendered with a dark matter profile centered at the barycenter of the three Local Group galaxies.

The number of particles assigned to each dark matter halo is determined as follows. A ref-

– 34 – erence number of dark matter particles, of order 105 or 106 depending on the desired resolution, is chosen for M31. The ratio of M31’s virial mass to this number gives the mass of each particle, which is constant across all objects for each run. The number of particles in the other objects is then proportional to their mass. The softening length must be adapted to the resolution. In general the mass resolution of dark matter particles must scale with the volume, which scales as 3 , where 1/3

 is the softening length. As a result, we vary  with the mass as  ∝ mparticle .

Combining the galaxies The four newly created objects are then combined by specifying their initial positions, velocities, spin parameters and inclinations. The halo spin parameter is dimensionless and defined in the following manner (Dutton & van den Bosch 2012): Jvir /Mvir p λhalo = √ fc 2Rvir Vvir

(28)

Here Jvir /Mvir is the specific angular momentum of the dark matter, where Jvir is the total angular momentum and Mvir the virial mass; Rvir is the halo’s virial radius, Vvir the circular velocity at the virial radius, and fc measures the deviation of the halo’s energy from that of a singular isothermal sphere truncated at the virial radius. fc is commonly set to 1 so that the spin parameter only depends on the halo’s mass and total angular momentum (Dutton & van den Bosch 2012). In our simulation, the Local Group halo, M31 and the Milky Way are not given any stellar matter, and as pure dark matter haloes they are spherically symmetric. In accordance with the method used for constructing the semi-analytic model, we ignore the inclination parameter (see Section 2) and arbitrarily set it to zero. Realistic portrayal of disk dynamics requires the presence of spin, and we assign all objects a spin parameter of 0.031, the average spin parameter of relaxed ΛCDM haloes being λhalo = 0.031 ± 0.001 (Dutton & van den Bosch 2012). In any case, the orbital sensitivity to the value of the spin parameter is very low. M31 is assigned stellar bulge and disk components according to the parameters listed in Table 1. Figure 7 shows face-on and edge-on views of the

– 35 – resulting synthetic stellar matter components for M33. The initial positions and velocities are the phase-space coordinates recorded at t = −5 Gyr in the matching semi-analytic run (see Figure 3).

Fig. 7.— Example of the synthetic stellar matter components generated with GADGET-2. Face-on and edge-on views of M33’s disk and bulge are shown with stars appearing color-coded yellow to red with increasing density.

Running the simulation The evolution of the self-gravitating particles in the four dark matter haloes and the stellar components is then simulated with GADGET-2 algorithms. In the low resolution trial runs, for each halo roughly 500,000 particles were tracked for a total of 5 Gyr, yielding a computation time of four hours on eight processors. Production runs were carried out in medium resolution, with dark matter particle masses of 1.6×107 M , yielding ∼ 105 dark matter particles each for the Milky Way and M31 haloes and accordingly lower/higher numbers for M33 and the Local Group halo (scaled according to their mass parameters µLG and µT ). Higher resolution is needed for the baryonic matter in order to better resolve potential disruptions of M33’s disk, so a stellar mass resolution of 2 × 105 M was chosen. A snapshot of each dark matter particle’s phase-space coordinates is recorded every 100 Myr.

– 36 – 4.

4.1.

RESULTS

Successful semi-analytic runs

µLG

µT

VW (km/s)

VN (km/s)

Run

0.0

1.0

-119

-21; -4.0

1

-98.5

-4.0

2

-37

-4.0

-119

-38; -21; -4.0

3

-98.5

-21; -4.0

4

-78

-4.0

-57.5

-4.0

-37

-4.0

-98.5

-4.0

0.25

0.5

1.0

1.15

5

Fig. 8.— Tree chart of the final group of successful parameter sets. µLG and µT are the Local Group and M33 mass parameters, vW and vN are the M31 proper motions in the western and northern direction. From the 375 trials in the second set of semi-analytic runs, only the 13 (µLG , µT , vW , vN ) combinations presented here yielded no collisions and successfully matched M33’s mass within the error bars. Five runs thought to be representative of different trajectory families were chosen and are highlighted in red in the VN column. In the remainder of this work they are labeled runs 1 through 5 from top to bottom of this chart.

– 37 – The parameter space outlined in Section 2.3 (see Table 3) was used for a first sequence of runs in search of successful parameter sets. No successful orbits were found for µLG ≥ 1 or µT ≥ 1.3, while solutions were found across the full range of transverse velocity components for M31. These findings warranted a second iteration of the semi-analytic model with a refined parameter space, as can be seen in Table 4. Orbits leading to M33 destruction due to close encounters with the Local Group halo were also eliminated in this iteration. This narrowed down the set of successful parameters to thirteen, which are shown as a tree diagram in Figure 8. For each of these arrays a reduced globular cluster population of twenty sample GCs was integrated and tracked over the course of 5 Gyr. In addition, a full M33 globular cluster population of 200 objects was simulated and projected on the plane of the sky, as described in Section 2.6. The globular clusters evolved as test particles in a background potential consisting of four dark matter haloes in addition to bulge and disk components from M31. Out of the thirteen possible parameter sets, we show five runs, each of them representative of a family of trajectories with similar behavior; these are shown in red in Figure 8. Three-dimensional images of the trajectories followed by each galaxy, as well as snapshots of the reduced globular cluster populations, can be seen in the upper left-hand panel of Figures 9 through 13. The trajectories of the three galaxies are shown in red for M31, blue for the Milky Way, and green for M33. The orbits are also projected onto the xy plane and the projections appear in paler colors. The center of the Local Group dark matter halo is shown as a black oval. Three snapshots of M33’s reduced globular cluster population are shown. The original distribution is plotted in dark green; the distribution at t = −2.5 Gyr in yellow; and the final, t = 0 distribution is in green. Figure 9 through 13 also show Hammer projections of the full GC populations overlaid on an infrared 2MASS map of the sky. Distinctive features of each population are better visualized in the zoomed-in insets in the corner of each sky map. M31’s and M33’s locations are represented by filled circles, in blue and green respectively; the circle size is not representative of the physical sizes of the galaxies. To provide a sense of scale, the M33-M31 separation on the sky is approximately 15 degrees, and the angular extent of M31 is ∼ 4◦ .

– 38 – 4.2.

Comparison of semi-analytic and GADGET results

Initial conditions recorded from the semi-analytic model at t = −5 Gyr were fed into the synthetic GADGET galactic components described in Section 3.2. The simulations were then run forward until t = 0. Decompositions of the resulting configurations in the xy, xz and yz planes are shown in Figures 14 through 18. In the GADGET runs, dark matter particles are color-coded yellow to black with increasing density. Past trajectories are shown as lines with the same colorcoding as Figures 9 through 13: M31 in red, the Milky Way in blue, M33 in green and the Local Group halo in black. In GADGET tracks are computed by calculating the center of mass of all particles belonging to a galaxy and tracking it in time. In the semi-analytic runs, the orbital tracks are calculated by recording the position of the point-particle galaxies at each moment in time. In the semi-analytic runs, starting locations are shown a black plus signs, and end positions as black squares. Despite the simplified treatment used in the semi-analytic model, we see that the tracks in the GADGET and semi-analytic runs have very similar behaviors. Not surprisingly, the most serious discrepancies occur for Run 5, where the Local Group halo is assigned a mass equivalent to the mass of M31, thereby strongly affecting the orbits of all three galaxies. Due to its large virial radius and initial location between M31 and the Milky Way, the Local Group halo is contiguous to the other galaxies, and point-mass treatment is no longer sufficient. As a result, the GADGET Run 5 predicts a close passage between the Milky Way and the Local Group halo, an event whose inevitable disruptive consequences have not been observed in our Galaxy. On this basis we can discard Run 5 as an implausible orbital history. In general the match between the GADGET and semi-analytic results is good, implying that given an initial parameter set, the semi-analytic model faithfully represents the Local Group’s orbital history.

– 39 –

Table 4.

List of parameters and their values.

µLG

µT

vW (km/s)

vN (km/s)

0.0 0.25 0.5 0.75 1.0

1.00 1.15 1.30 -

-119 -98.5 -78.0 -57.5 -37.0

-73.0 -56.0 -38.0 -21.0 -4.0

Note. — Parameter grid for second set of 375 runs.

– 40 –

Run 1

2MASS Image Mosaic: IR Processing and Analysis Center/Caltech & University of Massachusetts

Fig. 9.— Three-dimensional orbits, simulated M33 stellar matter and GC distribution for Run 1. On the overlaid map the GC distribution is shown as white dots, and a zoomed-in inset where the clusters appear as black crosses is provided in the upper right-hand corner.

– 41 –

Run 2

2MASS Image Mosaic: IR Processing and Analysis Center/Caltech & University of Massachusetts

Fig. 10.— Three-dimensional orbits, simulated M33 stellar matter and GC distribution for Run 2. On the overlaid map the GC distribution is shown as white dots, and a zoomed-in inset where the clusters appear as black crosses is provided in the upper right-hand corner.

– 42 –

Run 3

2MASS Image Mosaic: IR Processing and Analysis Center/Caltech & University of Massachusetts

Fig. 11.— Three-dimensional orbits, simulated M33 stellar matter and GC distribution for Run 3. On the overlaid map the GC distribution is shown as white dots, and a zoomed-in inset where the clusters appear as black crosses is provided in the upper right-hand corner.

– 43 –

Run 4

2MASS Image Mosaic: IR Processing and Analysis Center/Caltech & University of Massachusetts

Fig. 12.— Three-dimensional orbits, simulated M33 stellar matter and GC distribution for Run 4. On the overlaid map the GC distribution is shown as white dots, and a zoomed-in inset where the clusters appear as black crosses is provided in the upper right-hand corner.

– 44 –

Run 5

2MASS Image Mosaic: IR Processing and Analysis Center/Caltech & University of Massachusetts

Fig. 13.— Three-dimensional orbits, simulated M33 stellar matter and GC distribution for Run 5. On the overlaid map the GC distribution is shown as white dots, and a zoomed-in inset where the clusters appear as black crosses is provided in the upper right-hand corner.

– 45 –

Fig. 14.— Run 1 simulated by GADGET (left-hand side) and the semi-analytic code (right-hand side). The time elapsed since the beginning of the simulation is indicated in the upper left side corner of each panel. Particles are color-coded in function of density. In both panels the trajectory of each galaxy is shown in red for M31, blue for the Milky Way and green for M33. In the right-hand side panels, starting locations are shown as black plus signs, and end positions as black squares.

– 46 –

Fig. 15.— Run 2 simulated by GADGET (left-hand side) and the semi-analytic code from Section 2 (right-hand side).

– 47 –

Fig. 16.— Run 3 simulated by GADGET (left-hand side) and the semi-analytic code from Section 2 (right-hand side).

– 48 –

Fig. 17.— Run 4 simulated by GADGET (left-hand side) and the semi-analytic code from Section 2 (right-hand side).

– 49 –

Fig. 18.— Run 5 simulated by GADGET (left-hand side) and the semi-analytic code from Section 2 (right-hand side).

– 50 – 5.

DISCUSSION

Looking at Table 4 and Figure 8, we see that no solutions were found for µLG > 0.5 or for µT > 1.15. This outcome allows us to estimate upper bounds on the intergalactic dark matter and on M33’s initial mass. Run 5 was output by the semi-analytic model as a successful trajectory; comparison with GADGET output however indicates that it is not a valid parameter set, as it leads to a destructive encounter between the Milky Way and the Local Group halo. As a result, no successful orbital histories are found for µLG ≥ 0.5, yielding a first upper limit for the amount of intergalactic dark matter in the Local Group. We can also place a constraint on the original mass of M33 and the amount of tidal stripping that has taken place over the past 5 Gyr. The target mass for M33 at the end of each run is (2.2 ± 0.1) × 1011 M . Run 5 having been discarded by GADGET, no successful orbits were found for µT ≥ 1.15, that is, for MT,initial ≥ 2.53 × 1011 M . As a result, we estimate that a maximum ∼ 3 × 1010 solar masses of matter were stripped off of M33 over the past 5 Gyr. Figure 8 warrants a more detailed examination of the successful proper motion combinations. After eliminating Run 5, the remaining number of successful parameter sets is twelve. All solutions are for µT = 1.0. The allowed M31 proper motion range derived by Van der Marel & Guhathakurta (2008) is vW = −78 ± 41 km/s and vN = −38 ± 34 km/s. Importantly, looking at the rightmost column of Figure 8, we see that all runs with vN < −38 km/s were unsuccessful; that is, solutions are only found for the least negative values of proper motion in the northern direction. In fact, only one solution is found for the proper motion value at the center of the vN allowed range, corresponding to Run 3. Two-thirds of the twelve trajectories have vN = −4 km/s, which corresponds to the edge of the allowed proper motion range. These results indicate that, independently of the quantity of intergalactic dark matter, it is the less negative half of the allowed M31 vN range that produces plausible past orbits. Two-thirds of the successful sets are for µLG = 0.25, shown on the middle branch of Figure 8. These solutions span the entire range of proper motion in the eastern direction. Here again, a larger number of solutions is found for the more extreme values in the allowed range: vW = −119 km/s and vW = −98.5 km/s. For the less negative vW values, only the extreme value

– 51 – of vN = −4 km/s yields a solution. Similarly, in the top branch of the tree, the largest number of solutions is found for vW = −119 km/s. This indicates that M31’s velocity components are most likely to be found in the vN ≥ −38 km/s and vW ≤ −78 km/s ranges. We compare these new constraints on the range of allowed proper motion for M31 with Loeb et al. (2005), who derive their results by requiring the survival of M33’s stellar disk to tidal disruption from M31. Their simulations favor a proper motion amplitude of 100 ± 20 km/s and rule out the quadrant of positive vN and vW . The likely proper motion ranges derived here are consistent with this result. The 13 successful trajectory families shown in Figure 8 lead to an average proper motion amplitude of 93 ± 31 km/s, in agreement with the amplitude in Loeb et al. (2005). The timing argument discussed in Section 1 holds that the Milky Way-M31 orbit should be close to radial. Accordingly, we tested orbits with zero M31 transverse velocity with the semi-analytic model. The range of Local Group halo masses was from 0 to 1.2 Milky Way masses, and the range of initial M33 masses was from 1 to 1.8 times the current observed mass, 2.2×1011 M . Only one family of successful orbits was found. These all had initial M33 mass of 1.2 the current value, and Local Group halo masses ranging from 0.21 to 0.83 Milky Way masses. This implies that if the Milky WayM31 orbit is indeed radial, 0.2 to 0.8 Milky Way masses of intergalactic dark matter are required to explain the observed orbits, and approximately 20% of M33’s original mass must have been stripped due to interaction with M31. Darling (2011)’s discovery of water maser sources in M31 will lead to high-precision proper motion measurements in the coming years, effectively removing two dimensions of the parameter space in this work. These advances will immensely improve our understanding of past Local Group dynamics. We have found a general conservative upper bound for the intergalactic dark matter mass of approximately a Milky Way mass. Pinpointing the proper motion components of M31 will reduce the space of possible past galactic trajectories under different underlying dark matter distribution assumptions, thereby narrowing the range of possible dark matter masses.

– 52 –

Fig. 19.— Left panel: Stellar density map of the M31-M33 region from PAndAS (McConnachie et al. 2009). Dashed circles represent the maximum projected radii of 150 and 50 kpc from M31 and M33, respectively. Scale images of the disks of M31 and M33 are overlaid. Right panel: Map of globular clusters (black plus signs) resulting from Run 3, in equatorial coordinates. The green and blue dots represent the locations of M33 and M31 respectively.

We now move on to discussing possible signatures of the twelve orbital solutions on M33’s stellar matter. Figures 9 through 13 present snapshots of the final simulated stellar component of M33 in three different planes. In all runs, the original full symmetry is not preserved, and the stellar disk shows subtle density variations. This is visible for example in Figure 12 for Run 4, where in the xy plane the region towards negative x is overdense. For Run 3 (Figure 11) in the yz plane the disk appears warped. As another example, for Run 2 (Figure 10), the positive x side of the disk in the xz plane has been stretched. In this visualization the stellar disks however do not show strong structure akin to the features observed by McConnachie et al. (2009) in Figure 19. The north-south stream across M33’s center and the prominent stream in the south-west of M31 are Pan-M31 Archaeological Survey (PAndAS) discoveries. By mapping a vast panorama of the M31-M33 region of the sky, this survey seeks to confirm fundamental tenets of hierarchical galaxy formation. The globular cluster record can serve as a more useful discriminant for past orbits. Figure 20 shows maps in equatorial coordinates of final globular cluster population resulting from our semi-analytic model. The angular size of the globular cluster cloud is very large, extending much

– 53 – beyond the angular size of M33’s stellar disk, which is ∼ 1.2◦ (NED catalogue) (compare with the M31-M33 separation, which is ∼ 15◦ ). Prominent asymmetries are immediately obvious. Runs 1 and 4 yield distortion in the northeast - southwest direction. Run 5 preserves general symmetry but leads to accretion of 5 − 10 clusters onto M31. Run 2 shows strong east-west streams and additionally features a small group of very distant, isolated GCs in the southwest. The aspect of these features when projected onto the sky is shown in Figures 9 through 13.

Run 1

Run 2

Run 4

Run 5

Fig. 20.— Map of globular clusters (black plus signs) resulting from semi-analytic runs 1, 2, 4 and 5, in equatorial coordinates. The green and blue dots represent the locations of M33 and M31 respectively. In these plots north is up and east is right.

Run 3 (Figure 19) is the closest match to the available stellar data summarized in Figure 21. Tidal features more compact than in Runs 1 and 4 appear in the northeast - southwest direction, similar in scale to the stellar structure identified by McConnachie et al. (2009). In addition, Run 5 having been eliminated, Run 3 is the only solution presented here that yields accretion of GCs

– 54 –

Fig. 21.— Left panel: CFHT/MegaCam frames around M33 from Cockroft et al. (2011). Colored dots and boxes represent high-confidence globular cluster locations. The two circles represent projected radii of 10 kpc and 50 kpc centered on M33. Right panel: PAndAS map of the major substructures in M31 (Mackey et al. 2010). The solid black lines delineate survey coverage. Globular clusters are displayed as magenta points, and the areas bound in green indicate globular cluster overdensities.

onto M31. The ∼ 5 globular clusters predicted to the west of M31 coincide with the NW Group and the SW cloud in Figure 21. This image shows the vast extent in radial separation of M31’s outer halo clusters to the center of the galaxy. Mackey et al. (2010) make the argument that the alignment of many outer halo clusters with tidal streams represents a physical association and that, as a result, the majority of the remote globular cluster system has been assembled as a consequence of the accretion of cluster-bearing satellite galaxies. The investigation of the photometric, metallic and structural properties of M31 clusters by Huxor et al. (2011) provides further support for a scenario in which a significant number of M31’s outer halo GCs originate from accretion. Our results show that besides the accreted dwarf satellites present in M31’s halo, M33 could also be a source of accreted M31 outer halo clusters. Establishing the provenance of some of these globulars as coming from M33, as discussed in Section 2.5, would then provide strong evidence in favor of an

– 55 – orbit of type Run 3. The distinct features that appear in M33’s globular cluster cloud should allow distinguishing between various orbital solutions in order to identify the effects of Local Group dark matter and close M31 encounters. One caveat however is the large scale of the predicted globular cluster cloud as compared to the low number of observed outer halo clusters in M33. Cockroft et al. (2011) used the PAndAS M33 frames in search of outer halo star clusters, but found only one new GC in addition to the five previously detected by Huxor et al. (2009). This newly found cluster is shown in a double box in Figure 21, while the Huxor et al. (2009) detections are enclosed in a single box. The vast scale of this search, out to a projected radius of 50 kpc, is evident; yet only one new unambiguous globular cluster was confirmed. The authors suggest that this low number of M33 outer halo clusters compared to the 67 known outer halo GCs for M31 may be due to tidal stripping in a previous dynamical interaction with M31. The majority of the GCs predicted by the semi-analytic model lie within the frames of Cockroft et al. (2011)’s search and should presumably have been detected. Some of Cockroft et al. (2011)’s highest-confidence objects are awaiting confirmation from IR data in an upcoming paper; these may bring up the number of M33 outer halo clusters by at least a few. A number of steps could be taken to improve the dark matter mass constraints and stellar matter tracers presented in this work. On the semi-analytic side, including tidal stripping of M33 due to Local Group dark matter would be a more realistic constraint than the simple no-close-encounters condition currently in use. Further refining the grid spacing in a third and fourth iteration, notably in µLG space, would help improve the derived constraints on intergalactic dark matter, M33’s original mass, and M31’s proper motion. In our GADGET simulations, while a possibly interesting development, no runs have been carried out with M31 stellar matter because of the high resulting computational time. Similarly, it would be interesting and in theory possible to predict features in the globular cluster record with GADGET, as was done with the semi-analytic model. This would mean improving the mass resolution by a factor of 10 or more, such that dark matter particles would reach masses of ∼ 105 - 106 M , matching the approximate mass of a globular cluster containing a

– 56 – typical number of 106 stars. The dark matter particles in the M33 and M31 haloes could then be selected randomly, marked as globular clusters and tracked in time. Due to the presence of four dark matter haloes and the large scale of this simulation, this modification however would require large amounts of computational time. Finally, an interesting future line of investigation would be to look for possible quasar activity in M31 due to infalling M33 debris for orbits with particularly close past encounters between M33 and M31. The possibility of such AGN activity in the Local Group and consequences on galactic structure have not yet been considered. In the long-term this could be investigated by adding black hole and gas components to our GADGET simulations.

Many thanks to Avi Loeb and Laura Blecha for invaluable guidance, and to Mark Reid and Bob Kirshner for useful comments.

– 57 – REFERENCES Barkana R., Loeb A., 2001, PhR, 349, 125 Barmby P., Huchra, J. P., 2001, AJ, 122, 2458 Benjamin, R.A., 2008, in Massive Star Formation: Observations Confront Theory, eds. H. Beuther, H. Linz & Th. Henning, ASP Conference Series, Vol. 387, p.375 Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton Univ. Press, Princeton, NJ. Brunthaler A., Reid M. J., Falcke H., Greenhill L.J., Henkel C., 2005, Sci, 307, 1440 Brunthaler A. et al., 2011, Astron. Nachr., 332, 5, 461 Churchwell E., GLIMPSE Team, 2005, RevMexAA (Serie de Conferencias), 23, 53 Ciardullo R. et al., 2004, ApJ, 614, 167 Cockroft R. et al., 2011, ApJ, 730, 112 Corbelli E., Schneider S. E., 1997, ApJ, 479, 244 Courteau S., van den Bergh S., 1999, ApJ, 118, 337 Cox T. J., Loeb A., 2008, MNRAS, 386, 461 Darling J., 2011, ApJL, 732, L2 Dutton A. A., van den Bosch F. C., 2012, MNRAS, 421, 608 Eddington A. S., 1916, MNRAS, 76, 572 Forbes D. A, Bridges T., 2010, MNRAS, 404, 1203 Harris W. E., 1996, AJ, 112, 4. Catalog last updated in May-June 2011, accessed online on 12/08/2011 at http://spider.seds.org/spider/MWGC/mwgc.html Hernquist L., 1990, ApJ, 356, 359

– 58 – Huxor A. et al., 2009, ApJ, 698, L77 Huxor A et al., 2011, MNRAS, 414, 770 Johnson D. R. H., Soderblom D. R., 1987, /aj, 93, 864 Kahn F. D., Woltjer L., 1959, ApJ, 130, 705 Klypin A., Kravtsov A. V., Valenzuela O., 1999, ApJ, 522, 82 Kunkel W. E., Demers S., Irwin M. J., 1996, BAAS, 28, 931 Loeb A., Reid M. J., Brunthaler A., Falcke H., 2005, ApJ, 633, 894 Mackey A.D., Gilmore G. F., 2004, MNRAS, 355, 504 Mackey A. D. et al., 2010, ApJL, 717, L11 McConnachie A. W. et al., 2004, MNRAS, 350, 243 McConnachie A. W. et al., 2005, MNRAS, 356, 979 McConnachie A. W., et al., 2009, Nature Letters, 46, 66 McConnachie A. W., et al., 2010, ApJ, 723:1038. McLaughlin D. E., 1999, ApJ, 512, L9 Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493 Perrett K. M. et al., 2002, AJ, 123, 2490 Reid M. J. et al., 2009, ApJ, 700, 137 Rich R.M. et al., 2005, AJ, 129, 2670 Rohlfs K., Kreitschmann J., 1988, A & A, 201, 51 Sarajedini A., Mancone C. L., 2007, AJ, 134, 447

– 59 – Seigar M. S., 2011, preprint (arXiv:1103.3200) Shattow G., Loeb A., 2009, MNRAS, 392, L21 Sherwin B. D., Loeb A., O’Leary R. M., 2008, MNRAS, 386, 1179 Springel V., Yoshida N., White S. D. M., 2001, New Astronomy, 6, 79 Springel V., 2005, MNRAS, 364, 1105 Springel V. et al., 2005, Nature, 435, 629 Van der Marel R. P., Guhathakurta P., 2008, ApJ, 678, 187 Watkins L. L., Evans N. W., An J. H., 2010, MNRAS, 406, 264 Xue X. X. et al., 2008, ApJ, 684, 1143 Yoon S-J., Lee Y-W., 2002, Science, 297, 578

This preprint was prepared with the AAS LATEX macros v5.2.