Sculpting Andromeda – madetomeasure models for M31’s bar and composite bulge: dynamics, stellar and dark matter mass
Abstract
The Andromeda galaxy (M31) contains a box/peanut bulge (BPB) entangled with a classical bulge (CB) requiring a triaxial modelling to determine the dynamics, stellar and dark matter mass. We construct madetomeasure models fitting new VIRUSW IFU bulge stellar kinematic observations, the IRAC3.6 photometry, and the disc’s HI rotation curve. We explore the parameter space for the 3.6 masstolight ratio , the bar pattern speed (), and the dark matter mass in the composite bulge () within . Considering Einasto dark matter profiles, we find the best models for , and . These models have a dynamical bulge mass of including a stellar mass of (73%), of which the CB has (28%) and the BPB (45%). We also explore models with NFW haloes finding that, while the Einasto models better fit the stellar kinematics, the obtained parameters agree within the errors. The values agree with adiabatically contracted cosmological NFW haloes with M31’s virial mass and radius. The best model has two bulge components with completely different kinematics that only together successfully reproduce the observations (, , , ). The modelling includes dust absorption which reproduces the observed kinematic asymmetries. Our results provide new constraints for the early formation of M31 given the lower mass found for the classical bulge and the shallow dark matter profile, as well as the secular evolution of M31 implied by the bar and its resonant interactions with the classical bulge, stellar halo and disc.
keywords:
galaxies: bulges – galaxies: individual: Andromeda, M31, NGC224 – galaxies: kinematics and dynamics – Local Group – galaxies: spiral – galaxies: structure.1 Introduction
The Andromeda galaxy (M31, NGC224) is the closest neighbouring massive spiral galaxy, presenting us a unique opportunity to study in depth the dynamics of disc galaxy substructures, such as classical bulges and bars, the latter found in approximately 70 per cent of the disc galaxies in the local Universe (MenendezDelmestre et al., 2007; Erwin, 2017). In addition, our external perspective more easily proves a global view of M31 in comparison to the Milky Way, while as a similar mass disk galaxy, it allows us to place our home galaxy in context.
Historically, M31’s triaxial bulge has been mostly addressed as a classical bulge, while generally the bar component has been only qualitatively considered in the modelling of its stellar dynamics. However, an accurate dynamical estimation of the mass distribution of the stellar and the dark matter in the bulge must take into account the barred nature of M31’s central regions (Lindblad, 1956). More recent observations better quantify the triaxiality of the bulge which is produced by its box/peanut bulge () component (Beaton et al., 2007; Opitsch et al., 2017), a situation similar in many aspects to the Milky Way’s box/peanut bulge (Shen et al., 2010; Wegg & Gerhard, 2013; BlandHawthorn & Gerhard, 2016). The M31 is in addition entangled with a classical bulge () component (Athanassoula & Beaton, 2006). The is much more concentrated than the , with the two components contributing with 1/3 and 2/3 of the total stellar mass of the bulge respectively, as shown by Blaña et al. (2017, hereafter B17).
Each substructure in M31 can potentially teach us about the different mechanisms involved in the formation and the evolution of the whole galaxy. In particular, the properties of the component of M31 can give us information about the early formation epoch. Current galaxy formation theories consider classical bulges as remnants of a very early formation process, such as a protogalactic collapse, and/or as remnants of mergers of galaxies that occurred during the first gigayears of violent hierarchical formation (Toomre, 1977; Naab & Burkert, 2003; Bournaud et al., 2005). On the other hand, the massive of M31 provides us information about the evolution of the disc, as box/peanut bulges are formed later from the disc material. Box/peanut bulges in Nbody models are triaxial structures formed through the buckling instability of the bar, which typically lasts for , generating a vertically thick structure (Combes et al., 1990; Raha et al., 1991). Recent observations of two barred galaxies also show evidence of their bars in the buckling process (Erwin & Debattista, 2016). Box/peanut bulges are frequent being found in 79 per cent of massive barred local galaxies (, Erwin & Debattista, 2017). Note that box/peanut bulges are sometimes referred as box/peanut pseudobulges, however not to be confused with discy pseudobulges, which are formed by gas accreted into the centres of disc galaxies (Kormendy, 2013).
Moreover, on even longer timescales, box/peanut bulges and bars can interact through resonances with the disc and thereby redistribute its material, generating for example surface brightness breaks, as well as ringlike substructures (Buta & Crocker, 1991; Debattista et al., 2006; Erwin et al., 2008; Buta, 2017). Bars also transfer their angular momentum to the spheroid components, such as classical bulges (Saha et al., 2012, 2016), stellar haloes (PerezVillegas et al., 2017) and dark matter haloes (Athanassoula & Misiriotis, 2002), changing their dynamical properties. Furthermore, Erwin & Debattista (2016) show also with observations that classical bulges can coexist with discy pseudobulges and box/peanut bulges building composite bulges, a scenario that has also been reproduced in galaxy formation simulations (Athanassoula et al., 2016). This makes M31 a convenient laboratory to test formation theories of composite bulges and to better understand their dynamics.
To understand the formation and the evolution of Andromeda, and to accurately compare it with galaxy formation simulations, it is imperative to first determine the contribution and the properties of each of the substructures, such as their masses and sizes, as well as the dark matter distribution. In the outer disc region the gas kinematics constrain the dark matter distribution (Chemin et al., 2009; Corbelli et al., 2010). However, in the centre, the gas may not be in equilibrium due to the triaxial potential generated by the bar. Therefore, we model the stellar kinematics taking into consideration the triaxial structure of the . Opitsch (2016, hereafter O16) and Opitsch et al. (2017, hereafter O17) obtained kinematic observations of exquisite detail using the integral field unit (IFU) VIRUSW (Fabricius et al., 2012), completely covering the classical bulge, the and most of the projected thin or planar bar. In this paper we use these kinematic observations to fit a series of madetomeasure models that allow us to find constraints for the stellar and dark matter mass within the bulge region, as well as other dynamical parameters such as the pattern speed of the and the thin bar.
This paper is ordered as follows: Section 2 describes the observational data, its implementation, and the madetomeasure modelling of M31. Section 3 shows the results of the models that are separated in two main parts. In the first, Section 3.1, we present the main results of the parameter search exploration. In the second part, in Section 3.2, we present the properties of the best model and we compare it with the M31 observations. In Section 4 we conclude with a summary and a discussion of the implications of our findings.
2 Modelling the bulge of M31
Most dynamical models for the bulge of M31 assume a spherical or an oblate geometry for the bulge (Ruiz, 1976; Kent, 1989; Widrow et al., 2003; Widrow & Dubinski, 2005; Block et al., 2006; Hammer et al., 2010), making the mass estimations in the centre less accurate due to the barred nature of this galaxy. Nbody barred galaxy models can represent the bulge and the bar of M31 much better. However, finding an Nbody model that exactly reproduces all the properties of the M31 substructures is very difficult, because Nbody models depend on their initial conditions and on the bar formation and buckling instabilities, evolving with some degree of stochasticity. Therefore, here we use the Madetomeasure (M2M) method to model the bulge of M31 (Syer & Tremaine, 1996, hereafter ST96). This method can model triaxial systems and therefore it is the most suitable approach to model M31’s bar.
In the following sections we describe our technique that implements the M2M method to fit the kinematic and the photometric observations, which allows us to determine the main dynamical properties of the M31 composite bulge: the pattern speed of the bar (), the stellar masstolight ratio of the bulge in the 3.6 band () and the dark matter mass within the bulge ().
2.1 Madetomeasure method
We use the program nmagic that implements the M2M method to fit Nbody models to observations (De Lorenzi et al., 2007; De Lorenzi et al., 2008; Morganti & Gerhard, 2012; Portail et al., 2015; Portail et al., 2017a). In the original implementation of the M2M method (ST96) the potential and the model observables are calculated from the initial mass distribution of the particles, where their masses are then optimised to match observations, requiring a mass distribution of the particles that is close to the final model. In the nmagic implementation the potential is periodically recomputed to generate a system that is gravitationally selfconsistent.
A discrete model observable is defined for a system with particles with phasespace time () depending coordinates as:
(1) 
where is a known kernel that is used to calculate the distribution moments, is the weight of each particle that contributes to the observable, corresponding here to the particle’s mass. We increase the effective number of particles implementing an exponential temporal smoothing with timescale , obtaining the smoothed observable .
The observational data is composed by observations (e.g. number of pixels in an image), and by different sets of observations; here we work with one set of photometric observations and four sets of kinematic observations. Therefore, we generalise to observations with errors, and by observing the model similarly we have temporally smoothed model observables and kernels . The deviation between the model observables and the observations is defined by the delta
(2) 
and therefore the sum in time of is the chisquare of the temporal smoothed model observables and the observations.
The heart of the M2M method is the algorithm that determines how the weights of the particles change in time during the iterative fit to the observations. Here we use the “forceofchange” (FOC) defined by ST96 as:
(3) 
where is a constant adjusting the strength of the FOC. This relation is a gradient ascent algorithm that maximises in the space of the weights, defined in NMAGIC as
(4) 
Here the first term is just the total chisquare
(5) 
where , and are constants that balance the contributions between different sets of observables (Long & Mao, 2010; Portail et al., 2015). The term is an “entropy” introduced by ST96 that forces the weights of the particle distribution to remain close to their initial distribution, defined here as in Morganti et al. (2013); Portail et al. (2017a).
(6) 
where the “priors” are the averages of the weights of each of the stellar particle types. The entropy term also forces the model to slowly change its initial 3D mass density distribution. The factor balances the contribution between the entropy term and the chisquare term (De Lorenzi et al., 2007). Introducing the previous terms in equation 3 we have now the FOC equation
(7) 
With the observables that we define later the differential term becomes zero ().
2.2 Inputs to the M2M modelling from B17: initial Nbody model and projection angles
The M2M modelling requires an initial input particle model that contains the orbits required to construct a new model that successfully matches the observations. Therefore, we use the best matching particle model for the M31 bulge from B17, i.e. Model 1, which comes from a set of 72 Nbody models built with a box/peanut bulge () component and a classical bulge () component with different masses and scale lengths. These models evolved from a Hernquist density profile for the classical bulge and another for the dark matter halo, where none of these components have initial rotation. During these simulations the initial disc forms a bar that later buckles forming a , but leaving bar material in the plane which is the thin bar. The thin bar is aligned with the extending beyond this. We reserve the term “bar” for whole structure of the thin bar and the together. The bar and disc particles have the same label, as the bar evolved from the initial disc. The bar is entangled with the , where both structures evolve due to the transfer of angular momentum from the bar to the and the dark matter halo as well, gaining both rotation. The light of the bulge and the dominate in the centre, and therefore no stellar halo component is included. The number of particles used for the , bar and disc and the dark matter halo are , and .
Model 1 (see B17) has a concentrated with a 3D halfmass radius and a with a 3D semimajor axis of and a halfmass radius of . Within the radius , B17 measure a stellar mass of the composite bulge of , where the and the have and of the bulge total stellar mass, respectively. They estimate a stellar masstolight ratio in the 3.6 of . The initial dark matter halo mass within 51 is and within is . This model has a bar pattern speed of .
We tested our final results using another model from B17 as the input Nbody model for the M2M fits. This model had the same initial conditions as Model 1, except for the classical bulge mass being 30 per cent higher. We found only small differences in the final fitted M2M model.
We also need to project the M2M models on the sky to calculate the model observables defined later in Section 2.3, requiring the distance to M31 , the disc inclination angle , the disc major axis position angle , and the bar angle . For this we use the same quantities adopted as in B17: (McConnachie et al., 2005) (at this distance , and on the sky), (Corbelli et al., 2010), (de Vaucouleurs, 1958), and the bar angle measured in B17. The bar angle is defined in the plane of the disc (where the bar major axis would be aligned with the disc projected major axis for , see B17 Figure 1). Projecting into the sky results in an angle of measured from the line of nodes of the disc major axis, corresponding to a position angle of . We corroborate later in Section 3.1.5 that is the bar angle that best matches the photometry of the bulge, reproducing the bulge isophotal twist.
2.3 Fitting the photometry and IFU kinematics
In this section we describe how we prepare M31’s photometric and kinematic observational data to use as constraints for the M2M fitting with nmagic. The photometric data consist of an image of M31 from the Infrared Array Camera 1 (IRAC 1) . The kinematic data correspond to IFU observations of the bulge region of M31, and to HI rotation curves in the disc region. Consistently with the observations, we build model observables that measure the same quantities in the model and are used to fit to the equivalent data values. However, as we explain later in Section 2.8, to find our range of the best matching models we select a subsample of the fitted observations to compare them with the models. All the model observables defined here are temporally smoothed to .
2.3.1 Photometry I: IRAC 3.6 observations
The imaging data that we use come from the largescale IRAC mosaic images of M31 of the Spitzer Space Telescope (Barmby et al., 2006) kindly made available to us by Pauline Barmby. We use the IRAC 1 band that at 3.6 wavelength for two reasons: i) it traces well the old stars (bulk of the population) where the light is dominated by giant stars that populate the red giant branch (RGB), and ii) this band has the advantage of being only weakly affected by the dust emission or absorption (Meidt et al., 2014). The IRAC1 mosaic of Barmby et al. (2006) has pixels with size of 0.863 and covers a region of . We are interested in covering the inner bulge region, both the region where the dominates within in the projected radius, also where the is at in projection. We use a resolution of 8.63 () per pixel for the image, which is a convenient scale that faithfully shows the light gradients in the central region where the transition between the and the is. As we are interested in the scenario where the 10ring could be connected to the outer Lindblad resonance, we include the region of the stellar disc out to . We define an ellipse with this projected semimajor axis by fitting to the isophotes with the ellipse task in iraf. We mask the pixels of the image that are outside this 15 ellipse and proceed to fit the image. We also mask hot pixels in the image, foreground stars, and the dwarf galaxy M32. At the end of the filtering, the total number of photometric observable (pixels) used for the M2M fit is 170651.
The original image pixel values are in intensity . The surfacebrightness figures in the paper that are in are in the Vega system, and they are transformed from the original units using the 3.6 zeropoint calibration 280.9 (Reach et al., 2005). The conversion between the SB in and the luminosity is done using the absolute solar magnitude value (Oh et al., 2008), and multiplying by the pixel area .
We also require photometric error maps for the M2M modelling. Given that the M2M models are a representation of M31 in dynamical equilibrium, they cannot reproduce the observed substructures in M31 that are produced by perturbations such as spiral arms. Therefore, we include these smaller scale deviations between M31 and the models in the errors. For this we combined three types of error maps: the observational error , the variability between pixels and the asymmetry error . The first error () is calculated from the square root of the sum in quadrature of the pixel error and the standard deviation for each pixel that comes from the original 0.863 pixels. The typical errors are between one and 5 per cent of the intensity depending on the pixel location in the image. This error is smaller than the variability observed between contiguous pixels and so we therefore include a second error that takes into account the pixeltopixel scatter. The surfacebrightness image of our M2M models is smoother than the observations. We take into account this variability by including in the photometric error the standard deviation within a radius of one 8.63pixel around each pixel of the image, obtaining the error . Finally we also include the variability observed at kiloparsec scales due to substructures like the spiral arms beyond the bar region, and the 10ring. For this we subtract the image with the same image, but rotated 180°around the centre of the bulge, obtaining . The bulge is roughly symmetric making this term smaller in the bulge than in the disc region. The combined photometric error per pixel with is then:
(8) 
2.3.2 Photometry II: model observables and the masstolight ratio ()
The photometric model observables consist of an array of pixels that extends from the bulge centre out to the disc until 15 along the disc major axis, where each model pixel uniquely corresponds to each observed pixel, with the same pixel size (). Each th pixel measures the stellar masses of particles that pass through each pixel, which are converted to light in the 3.6 band using the stellar masstolight ratio . The total light per pixel is the photometric model observable with :
(9) 
where the light per particle () is just . We define three masstolight ratio parameters in the 3.6 band: for the classical bulge, for the and for the outer disc, which are assigned to the particles according to the relation:
(10) 
where the particles are assigned everywhere, and the bar and disc particles at the cylindrical radius are assigned within , and if they are outside this radius. The last Gaussian term provides a smooth transition of from the value of to the value in the disc , where is the transition radius between the end of the thin bar and the disc (B17), and is the scale of the transition ().
In Section 2.8 we explain in more detail the different masstolight values that we explored, however in our fiducial M2M fits we assumed . In Sections 3.1.2 and 3.1.3 we explore further different values for each component, finding only small differences compared with our range of best models. From equation 9 we have that the photometric kernel () is
(11) 
2.3.3 Kinematics I: M31 Bulge IFU observations
O16 and O17 obtained kinematic IFU observations of the central region of M31 using the McDonald Observatory’s 2.7meter Harlan J. Smith Telescope and the VIRUSW Spectrograph (Fabricius et al., 2012). They cover the whole bulge and bar region and also sample the disc out to one disc scale length along six different directions, obtaining lineofsight velocity distribution profiles (LOSVDs). From this they calculate the four GaussHermite expansion coefficient moments (Gerhard, 1993; Bender et al., 1994), and obtain kinematic maps for the velocity , the velocity dispersion and the kinematic moments and . The velocity maps are corrected for the systemic velocity of (de Vaucouleurs et al., 1991). Note that the light weighted mean lineofsight velocity and the light weighted velocity standard deviation (or dispersion) , differ slightly from and when the LOSVDs deviate from a Gaussian distribution ( or or nonzero higher moments). This is because and are instead chosen so that the lower order GaussHermite terms, and , are zero.
We regrid the kinematic observations into new maps with the same spatial resolution of the photometric data. The new values of , , and are calculated from the error weighted average of the original values, leaving 13400 measurements for each kinematic variable, and therefore 53600 kinematic values in total. The regridded observational kinematic errors (, with ) are calculated from the standard deviation of the error weighted average. Similarly to the photometry, we combined the new observational error and the error due to the variability between different kinematic pixels within one pixel radius (), obtaining a total kinematic error per observable and per set of:
(12) 
2.3.4 Kinematics II: model observables
We now proceed to build the kinematic model observables. Because the kinematic observations are performed in the V band, we need to include the effects of dust in our model observables. A further description is given later in Section 3.2.3. Our dust absorption implementation consists of using M31 dust mass maps (Draine et al., 2014) converted to a V band absorption map by the dust model of Draine & Li (2007)
(13) 
We convert this to a 3D absorption map , deprojected as
(14) 
where for simplicity we assume that the dust is located in the plane of the disk, and therefore
any stellar th particle that is temporarily passing behind the disc at the moment that the kinematic model
observable is measured, is attenuated by the corresponding value of in the th pixel.
So that the kernel of Equation 1 does not depend on weight, we desire kinematic model observables that are linear in the particle weights. Therefore, we fit the GaussHermite moments of the observations, and , instead of directly fitting and (De Lorenzi et al., 2007). The model kinematic observables are then the lightweighted GaussHermite coefficient moments, calculated as in De Lorenzi et al. (2007), but with the inclusion of dust absorption:
(15) 
Here , and are the dimensionless GaussHermite functions (Gerhard, 1993),
(16) 
where are the standard Hermite polynomials, are
(17) 
where is the particle’s lineofsight velocity. The expansion is performed with the observational values of and so that while 1 and 2 are zero in the observations, they are in general nonzero when observing the model. From this we obtain the light weighted model observables , , , and . The corresponding kinematic kernel that changes the weights of the particles is
(18) 
Concordantly, the observational data that we fit are the GaussHermite moments , , and , which are lightweighted by the extincted light model observable
(19) 
This is then used to light weight the kinematic observations e.g. , obtaining the observations that we fit: , , and .
The errors for and are calculated from the observations and as in van der Marel & Franx (1993); Rix et al. (1997).
(20) 
Then, the kinematic errors , , and are also lightweighted in the form
,
which gives larger errors to the regions with more light extinction.
From this we obtained the light weighted errors , , , .
We also test our best model fit considering no dust absorption () and
a constant value .
To facilitate sidebyside comparison of the model with the observations, and also for the selection of the range of best models defined in Section 2.8, we also compute after the M2M fitting the temporally smoothed and of the model, and use these values to calculate and of the model. For this we observe the model and calculate , , , of the model using equation 15, but in equation 17 we replace and of the observations by the mean velocity and the velocity standard deviation of the model. The nonlight weighted quantities are recovered dividing by , i.e. h1=H1/ and similarly for , , . The parametrisation of the LOSVD with the GaussHermite moments dictates that the variables and are chosen such that and are zero. If this is not the case we use again the approximation (van der Marel & Franx, 1993; Rix et al., 1997) to correct and replace the old values of the velocity and the dispersion () with the new values () that result in new h and h values closer to zero:
(21a)  
(21b) 
We repeat the previous corrections observing the model and calculating the new , , , from the new dispersion and velocity using equation 15, repeating this iteratively until the terms and converge to zero or values smaller than the observational errors.
2.4 Adjusting the dark matter mass within the bulge (), and fitting the HI rotation curve
Our goal is to determine the dark matter mass within 3.2 of the bulge , by exploring a vast range of values given in Section 2.8. For this we change the initial dark matter mass distribution of the input Nbody model to match a target analytical profile. As we also want to explore the cusped or cored nature of the dark matter density in the central region, we consider different shapes for the target dark halo, making M2M models with two different target profiles. We consider the Einasto density profile (Einasto, 1965) which has a central core, parametrised here as:
(22) 
where is the elliptical radius for a flattening , is the scale length, is the central density and is the steepness of the profile. We also comte models with a NavarroFrenkWhite (NFW) dark matter mass density profile, which has a cuspy central profile (Navarro et al., 1996), parametrised here as
(23) 
where is the central density and is the scale length.
The parameters of these target analytical profiles are determined during each M2M run similarly to Portail et al. (2017a), by fitting the dark matter halo profile together with the current stellar mass distribution to match: i) the dark matter mass enclosed within an ellipsoidal volume of the major axis of the bulge () is fixed to the chosen value , with (or from the particles); and ii) that the total circular velocity of the model matches well the disc HI rotation curve data (Corbelli et al., 2010) described in Section 2.4.1.
To adjust the particle dark matter distribution to the target analytical dark matter profile we also use the M2M method (De Lorenzi et al., 2007). This is done by expanding the initial dark matter density distribution of the particles and the target analytical dark matter density profile in spherical harmonics, which are then fitted with the M2M scheme. The adaptation of the dark matter particles is performed while the photometric and the stellar kinematic observations are also being fitted.
A change in the dark matter mass profile may significantly change the total circular velocity, particularly in the disc region, affecting the orbits of the particles. This is not desirable for particles in the disc that should remain on nearcircular or epicyclic orbits. To alleviate this we measure the circular velocity for a th particle before and after the potential update, and then rescale the velocity of the particle living in the old potential to a new velocity given by the new potential by multiplying its velocity by the factor that is the ratio between the new and the old circular velocities:
(24) 
using the spherical radius vector for the dark matter and particles that have a spheroidal geometric distribution, and the cylindrical radius for the disc particles.
2.4.1 Kinematics III: HI rotation curve
We use the deprojected azimuthally averaged HI rotation velocity curve estimated by Corbelli et al. (2010) to fit the total circular velocity of our M2M models modifying the dark matter profile for a given (see Section 2.4) . This data extend from 8.5 out to 50. We do not fit the rotation curve beyond 20, for two reasons: i) the contribution of the mass of the HI disc to the circular velocity beyond this radius becomes as important as the stellar disc (Chemin et al., 2009), and ii) the outer disc shows a warp () changing the inclination with respect to the inner part of the stellar and gaseous discs (Newton & Emerson, 1977; Henderson, 1979; Brinks, E.; Burton, 1984; Chemin et al., 2009). This region includes the 10ring and the 15 ring structures Gordon et al. (2006); Barmby et al. (2006). We do not include the mass of the gas component in the potential as the gas mass and surface mass contribution within 20 is estimated to be less than 10 per cent of the stellar mass () (Chemin et al., 2009) and, as we show later in Section 3.1.3, the choice of different dark matter profiles introduces variations larger than this.
2.5 Bar pattern speed adjustment ()
The pattern speed of the bar of the model found in B17 is . As we want to find constraints for this quantity, we also explore pattern speeds (see Section 2.8). To change the initial pattern speed, we adiabatically and linearly change its initial value to the desired final value with a certain frequency defined in Section 2.7 (see MartinezValpuesta, 2012; Portail et al., 2017a). This pattern speed change is performed while the kinematic and the photometric observables are fitted and the potential is frequently recalculated from the new density distribution, resulting at the end of the M2M fit in a selfconsistent dynamical system.
2.6 Potential solver and orbital integration
As in Portail et al. (2017a), the NMAGIC modelling here uses the hybrid particlemesh code from Sellwood et al. (2003) to calculate the potential from the particle mass distribution. The potential solver uses a cylindrical mesh Fourier method to calculate the potential for the disc and the bulge components (Sellwood & Valluri, 1997). Due to the disc geometry and our interest in resolving the vertical and the in plane distribution, instead of using a spherical softening, we use an oblate softening with 67 in the plane and 17 in the vertical direction. The potential of the particles of the dark matter component is calculated using a spherical mesh with a spherical harmonics potential solver that extends to 42 and includes terms up to the 16th order (De Lorenzi et al., 2007). The cylindrical mesh extends in the disc plane out to and in the vertical direction, and any stellar mass particle that extends beyond the limits of this mesh is considered during the run in the spherical mesh for the calculation of the potential.
The orbits of the particles are integrated forward in time with an adaptive leapfrog algorithm using the acceleration due to the gravitational potential of all the particles. In the nmagic M2M implementation the rotating bar is kept fix in the reference frame of the potential by rotating the phasespace coordinates of all the particles around the axis at the same rate of the pattern speed of the bar, but opposite in sign (MartinezValpuesta, 2012; Portail et al., 2015) (note that the rotated system is still in an inertial frame).
The integration time is measured in iteration units , with a time step of (see B17, ). We require that the orbits always have at least 1000 steps per orbit.
2.7 M2M fitting procedure and parameters
Each M2M fitting done here with nmagic takes a total number of iterations of , where each fit is divided in three main phases. The first phase uses , and is when the temporal smoothed measurements of the model observables are calculated. The temporal smoothing scale is , and it is chosen to be larger than the period () of a circular orbit at 5 with circular velocity , which typically is .
The second phase is when the M2M fitting is performed, and it takes . The bar pattern speed is adjusted during this phase, starting at and finishing at , with an update of the new value every . During the second phase the total mass of the system may change. Therefore, we recalculate and update the potential from the new mass density distribution every . These regular potential updates are important to build a system that is gravitationally selfconsistent with its density.
The final phase is the stability check that takes , where the M2M fitting stops and the model is only observed. During this phase we recover the values of , , and for the model according to equation 21 correcting them every .
The FOC parameters , and of equation 7 are chosen sequentially. We first fit only the photometry, leaving the parameters and fixed to zero and varying only (the parameter , or , normalises and for simplicity is set to ). We measure the reduced chisquare (terms in equation 5), for the photometry () in the bulge region finding the relation between and shown in Figure 1 in the top panel. For too small the photometry does not have the power to change the model and so the is large. For too large the photometry has too much power, changing the particle weights too quickly compared to the orbital timescale, so that only the local observable (or pixel) that the particle is crossing is fitted. An optimum value allows the weight to change an averaged amount once it crosses all the observables that are along the particle’s orbit, so that its weights converge to a constant value. We find this optimum value at the minimum , when .
Second, we find the best for the IFU kinematic observables (where ) defined also as . We use the previous best and fit the photometry together with the IFU kinematics for several values of . We measure both the photometric and kinematic reduced chisquare in the bulge region, obtaining the relations versus shown in Figure 1 (second panel). The photometric has smaller values for small , and would be minimized for , because then only the photometry would be fitted without the additional kinematic constraints. As increases, the kinematic observations have more power tailoring the model towards fitting the kinematics as well, as we see the kinematic decreasing for larger , which worsen the photometric if the kinematics get too much power. Similarly to , if increases too much, both the photometric and the kinematic get worse. We find an optimal value of for the minimum kinematic while the photometric is still small.
To find the best parameter for the dark matter halo fitting, we fix the previously found parameters and and test different values of versus the reduced chisquare of the dark matter halo density (Figure 1 third panel). We find the minimum at , where the photometric and the kinematic remain almost unchanged.
We determine the entropy magnitude term to be (Figure 1 bottom panel) in the same way, fixing the previous parameters and choosing the largest that still has small values for the photometry and the kinematics.
After setting the fitting parameters, we run M2M fits, showing an example in Figure 2 where the reduced chisquares of the model observables from equation 5 are plotted versus time (iterations). In the phase the model temporal smoothed observables are calculated decreasing at first and then staying constant. Then the fitting phase starts where of the photometry, kinematics and the dark matter halo decrease in time. Finally in the stability check phase the values of increase slightly.
2.8 Exploring the effective potential parameters
B17 find good constraints for the mass ratio between the and the while the dark matter distribution and the bar pattern speed are less constrained. Here, we use stellar kinematic and photometric observations as targets to better determine these properties. While the M2M method has the power to change the orbital distribution, thereby changing the model , , and to match the observed kinematics, there are macroscopic potential parameters that limit the orbital phase space, and therefore a particular model will fit the data as well as these macroscopic parameters allow. Here we have three important dynamical quantities that are inputs to the M2M modelling and that impact the effective potential: the pattern speed of the bar , the stellar masstolight ratio of the bulge in the 3.6 band which, for a well fitted target observed luminosity, determines the total stellar mass in the bulge , and the amount of dark matter in the bulge region . Therefore, we need to apply a method of metaoptimization where each M2M model is an optimisation itself that finds the orbit distribution that best matches the observations for fixed potential parameters. Then we vary , and around reasonable values that we estimated from the literature, and then we find the range of values for which the M2M models overall best reproduce all sets of observations. To explore these three global parameters we create one cube (or grid) of model parameters for the Einasto dark matter profile, and a second cube for the NFW dark matter halo profile, where each model has the coordinates:
(25) 
For the Einasto cube we explore in the range of in steps of to produce a low resolution grid that allows us to quickly find the best fitting region, and then we include more values between in steps of . For we explore in steps of . For we explore the range in steps of , building then a cube of parameters with , i.e., a total of 1040 M2M models with the Einasto dark matter profile.
For the NFW cube we explore in the range in steps of . For we explore in steps of . For we explore in steps of , giving a cube of parameters with , i.e., a total of 420 M2M models for the NFW cube.
Dark matter haloes are expected to be flattened in the central part of disc galaxies due to the influence of the disc gravitational potential. Widrow et al. (2003); Widrow & Dubinski (2005) explored different flattening values for the dark halo of M31, finding reasonable fits between and 1.0. Here we use a dark halo flattening of as our fiducial value for both dark matter density profiles, but we test the effects of different values on the final results. We explore and , finding stellar mass distributions for the disc and the central region of the similar to the fiducial model. This is discussed further in Section 3.1.3.
2.9 Selection of bestmatching models in effective potential parameter space
The selection of the bestmatching models from the parameter grid just discussed cannot be done by straightforward minimization, because with the extended, highquality data available here, systematic effects play a dominant role. These include uncertainties in the dust modelling, intrinsic asymmetries in the observed surface brightness distribution (Figure 17 in Section 3.2.2 below), uncertainties in the parametrisation of the dark matter density distribution, and likely gradients in especially between the BPB and adjacent disk. Because of these systematic effects no model is found to give the best fit simultaneously in all regions of M31, and both photometric and kinematic observables.
In addition, while the M2M models are fitted to an impressive number of 224251 photometric and kinematic data values (pixels), the spatial distributions of photometric and kinematic pixels and their residuals are substantially different. (i) Typical errors can differ between different variables, e.g., between and , or and , leading to different ranges of ; see Figure 3. (ii) For the same variable set, the errors depend on the spatial regions considered; e.g., relative photometric errors are smaller in the central bulge than in its outer parts or in the disc region (Figure 4 below). Yet in all of these locations the data may contain signatures important for specific physical properties of the system.
In consequence, combining all values linearly in one total and finding the M2M model with that minimum total chisquare will not adequately capture the entire structure of M31; e.g., it will lead to a model providing a good fit of the region, but to an unsatisfactory fit in the smaller central region. In the following we therefore describe an alternative procedure which we believe leads to a more robust selection of the overall bestmatching models for M31 given the available data.
2.9.1 Building a metric for the comparison with the observational data: five chisquare subsets
Separating the observables , , , we first build five subsets of normalized . These are motivated by the properties of the system that we are modelling, which is built from the three main substructures , and disc that we want to fit simultaneously well. The dominates the light in M31 within , defined as region CBR. Further out the dominates the light within ellipses with semimajor axis (region BPR), and even further out the disc dominates (Beaton et al., 2007, B17); see Figure 4. We define five subsets of normalized s:

central photometry (): we measure the normalized (per data point) of the photometry () in the inner CBR within a diameter of 40 (150, ). With this we search for models that match the cuspy light profile in the centre of M31’s bulge.

photometry (): we measure the normalized of the photometry in region BPR (Figure 4).

dispersion (): B17 show that the and the of M31 have different kinematic properties. Hence, we calculate the normalized of the dispersion only in the BPR. This allows us also to find the dynamical mass within the bulge.

mean velocity (): Tremaine & Weinberg (1984) showed that the bar pattern speed is related to the LOS velocity () and the photometry. Therefore, we constrain the bar pattern speed with the normalized of the mean LOS velocity in the bar region BPR.
In this way, each model is evaluated with five normalized parameters, . While the GaussHermite coefficients and and all observables in the disc region are also fitted in each of the M2M models, we do not include subsets for them in the best model selection; later we show that the best models selected by the five subsets defined above satisfactory reproduce these observations as well.