Dynamical Models for M31 - Bulge & Bar

Andromeda chained to the Box -- Dynamical Models for M31: Bulge & Bar


Andromeda is our nearest neighbouring disk galaxy and a prime target for detailed modelling of the evolutionary processes that shape galaxies. We analyse the nature of M31’s triaxial bulge with an extensive set of N-body models, which include Box/Peanut (B/P) bulges as well as initial classical bulges (ICBs). Comparing with IRAC 3.6 data, only one model matches simultaneously all the morphological properties of M31’s bulge, and requires an ICB and a B/P bulge with 1/3 and 2/3 of the total bulge mass respectively. We find that our pure B/P bulge models do not show concentrations high enough to match the Sérsic index () and the effective radius of M31’s bulge. Instead, the best model requires an ICB component with mass and three-dimensional half-mass radius . The B/P bulge component has a mass of and a half-mass radius of . The model’s B/P bulge extends to in the plane of the disk, as does M31’s bulge. In this composite bulge model, the ICB component explains the velocity dispersion drop observed in the centre within , while the B/P bulge component reproduces the observed rapid rotation and the kinematic twist of the observed zero velocity line. This model’s pattern speed is , placing corotation at . The outer Lindblad resonance (OLR) is then at , near the 10-ring of M31, suggesting that this structure may be related to the bar’s OLR. By comparison with an earlier snapshot, we estimate that M31’s thin bar extends to in the disk plane, and in projection extends to .

galaxies: spiral — galaxies: bulges — galaxies: bar — galaxies: kinematics and dynamics — galaxies: individual (Andromeda, M31, NGC224)) — galaxies: Local Group – galaxies: structure — methods: nbody

1 Introduction

The Andromeda galaxy (M31, NGC224) is the nearest spiral galaxy, located at a distance of from the Milky Way (MW) (McConnachie et al., 2005), which makes it one of the best cases to study in great detail the structural and evolutionary properties of spiral galaxies. It shows several signs of a hierarchical formation history through its satellites, the Giant Stream and its old stellar halo suggesting a history of minor mergers (Ibata et al., 2001; Tanaka et al., 2010) or also a possible major merger (Hammer et al., 2010; Bekki, 2010). It also shows signs of secular evolution such as a massive star forming disk, spiral arms and rings (Gordon et al., 2006; Barmby et al., 2007; Chemin et al., 2009). The massive bulge of Andromeda has commonly been considered in the literature as a classical bulge, or more recently, as a classical bulge with “pseudobulge trimmings” (Mould, 2013).

According to theory, classical bulges are considered as elliptical galaxies in the centres of disks, as remnants of a very early formation process, or as remnants of mergers of galaxies occurred during the first gigayears of a violent hierarchical formation (Toomre, 1977; Naab & Burkert, 2003; Bournaud et al., 2005; Athanassoula et al., 2016). Another proposed mechanism to form early bulges consists of instabilities in the early disk, forming clumps and star clusters that spiral in due to dynamical friction, merging later and forming a bulge (Noguchi, 1999; Immeli et al., 2004; Bournaud et al., 2007). The comprehensive review of Kormendy (2013) describes in general how disc galaxies can also develop bulges through secular evolution, referred also as pseudobulges. The term pseudobulge groups disky bulges and box/peanut bulges together. In this scenario pseudobulges are a manifestation of the evolution of the disk, contrary to the classical bulges.

Fisher & Drory (2008) use a sample of spiral galaxies and their photometric properties such as the Sérsic index , the effective radius , and the disk scale length , to distinguish pseudobulges from classical bulges. They conclude that usually pseudobulges have Sérsic indices and classical bulges indices . Fabricius et al. (2012) additionally show a correlation between the Sérsic index, the velocity dispersion and the rotation. They compare the Sérsic index and the velocity dispersion averaged within one tenth of the effective radius . Their sample shows that classical bulges, defined morphologically or through a Sérsic index of , tend to have higher , rarely getting lower than km s and with a mean of the sample of km s. Pseudobulges, defined with , show lower , some as low as km s, and a mean value of km s for the sample. However these classification criteria have only statistical meanings, as a particular galaxy, and its bulge, may actually present properties of both types of bulges.

Simulations of wet mergers by Keselman & Nusser (2012) showed that classical or primordial bulges may also present characteristics of pseudobulges, such as rotation and , which makes even harder to distinguish the origin of its properties. Furthermore, Erwin et al. (2014) showed with observations that classical bulges can coexist with disky pseudobulges and box/peanut bulges, a situation that has also been reproduced in simulations (Athanassoula et al., 2016).

Saha et al. (2012) showed that if a low mass non-rotating classical bulge is present during the formation and evolution of the bar it absorbs angular momentum from the bar (see also Hernquist & Weinberg, 1992; Athanassoula, 2003), and subsequently the bulge manifests cylindrical rotation like the bar, thereby making the final combined structure with the bar difficult to disentangle. However, more massive classical bulges do not reach cylindrical rotation (Saha et al., 2016).

Observations show that M31’s bulge is triaxial, with a photometric twist relative to the disk (Lindblad, 1956; Hodge & Kennicutt, 1982; Beaton et al., 2007), and the misalignment of the bulge kinematic major axis and its photometric major axis (Saglia et al., 2010). A known dynamical mechanism to generate a triaxial structure is the buckling instability of the bar, which generates the box/peanut (B/P) bulge in N-body models (Combes et al., 1990; Raha et al., 1991, see Section 3.1.2) . Erwin & Debattista (2016) have recently found indications that two barred galaxies (NGC 3227 and NGC 4569) are currently in the buckling phase.

Several hints such as: the kinematics of the bulge, the boxiness of the surface-brightness contours, and the twist of the photometric major axis of the bulge’s boxy region and the disk’s major axis, suggest that M31 could contain a B/P bulge and therefore may also host a thin bar. The inclination of (Corbelli et al., 2010) of the disk presents challenges, but the proximity of M31 and the excellent available photometric and kinematic data (Barmby et al., 2007; Chemin et al., 2009; Corbelli et al., 2010; Saglia et al., 2010), make it possible to explore this scenario. Courteau et al. (2011) (hereafter Co11) produced photometric models from the IRAC 3.6 data and find a bulge-to-total luminosity ratio of a Sérsic index of , an effective radius or (M31 distance implies a unit conversion of , and ), and a disk scale length of , giving a ratio of . These values would locate M31’s bulge just between the classical and the pseudobulges, according to Fisher & Drory (2008). However, the kinematics in the bulge region (Saglia et al., 2010, see the major axis velocity and dispersion profiles in their Fig.3), would classify it among the group of classical bulges (Fabricius et al., 2012), due to the high velocity dispersion at of km s and the rather slow rotation in the bulge region (km s at ) .

Another interesting hint of the bulge type of M31 comes from the metallicity analysis and the age estimation of the stars. Saglia et al. (2010) find the presence of a very old stellar population in the central part within , with an average age of 12, which would support the idea of a classical component in the bulge. Meanwhile, from outwards the mean age of the stellar population drops to 8 Gyr. This is compatible with the isophotal analysis performed by Beaton et al. (2007), where the authors present strong evidence for the presence and projected extension of the B/P bulge in M31, by quantifying the morphological properties of the isophotes of M31 2MASS 6X images in the J,H,K bands. Within the round isophotes resemble a classical bulge with low ellipticity. Beyond the boxy isophotes emerge, also increasing the ellipticity, extending until , where the strongest boxiness is at . Athanassoula & Beaton (2006, hereafter AB06) used four N-body models of barred galaxies based on Athanassoula (2003), two with a pure B/P bulge and two with a classical bulge combined with a B/P bulge, to compare qualitatively the shape of their isophotes with the isophotes of M31 in the J band, concluding the presence of a dominant B/P bulge and also a classical bulge component in M31.

With the IRAC 3.6 data that are less affected by dust absorption than the 2MASS 6X data and go deeper into the bulge and disk, together with the new kinematic data of Saglia et al. (2010); Chemin et al. (2009); Corbelli et al. (2010) and Opitsch et al. (in preparation), we can quantitatively compare M31 with -body models. From these observational quantities we can build an understanding of the inner structure of M31’s bulge and estimate the contribution of each possible subcomponent i.e. a classical bulge, a B/P bulge and a thin bar. The goal of this paper is to explore dynamical models for the bulge of M31 that consist of two components, a classical component and a B/P bulge, and also test pure B/P bulge models.

The paper is ordered as follows: Section 2 characterises the observational data, Section 3 describes the set-up of the simulations, and our technique to compare the simulations with the observations. Section 4 shows simultaneously the results for M31 and for the simulations. The morphological and photometric analysis are presented in Sections 4.2 and 4.3. In Section 4.5 we present the properties of the best model. In Section 5 we discuss properties of triaxial models in the literature, and in Section 6 we conclude with a summary of our findings and the implications for M31.

2 Observational data: M31 IRAC 3.6 image

The imaging data we use for our analysis come from the large-scale IRAC mosaic images of the Spitzer Space Telescope (Barmby et al., 2007) – specifically, the 3.6 IRAC1 mosaic shown in Fig.2, kindly made available to us by Pauline Barmby. The near side of the disk is in the upper part of the figure, pointing north-west, as evidenced by the dust lanes and reddening maps (Walterbos & Kennicutt, 1988). Careful inspection of the edges of the mosaic showed extended regions with slightly negative pixel values, suggesting a small oversubtraction of the sky background. Although this has negligible effect on our analysis of the bar and bulge region, we corrected the image by measuring the background in the regions furthest from the galaxy centre along the minor axis. The original IRAC1 mosaic of Barmby et al. (2007) has a pixel scale of 0.863 arcsec/pixel and a size of pixels ( degrees), using here for the analysis a resolution of 8.63/pixel. We fit ellipses to the 3.6 isophotes using the ellipse task in the stsdas package in iraf; this task is based on the algorithm of Jedrzejewski (1987). The algorithm fits ellipses to isophotes by minimizing the deviations in intensity around a given ellipse; for an ellipse with a given semi-major axis , this is done by iteratively adjusting the ellipse centre (pixel position in and ), orientation (position angle ), and ellipticity ( or , where is the semi-minor axis of the ellipse). From ellipse we obtained the azimuthally average (AZAV) intensity () profile. All magnitudes are in the Vega system, and we use the 3.6 absolute solar magnitude (Oh et al., 2008). We convert the intensity profile to surface-brightness profile (SB) using the 3.6 zero-point calibration Jy (Reach et al., 2005). For each fitted ellipse, the algorithm also expands the residual variations in intensity around the ellipse in a Fourier series:


where is the mean intensity along the best-fit ellipse; the and components are automatically zero for the best-fit ellipse. In fact, the actual computation divides the and higher coefficients by the local radial intensity gradient and the ellipse semi-major axis to generate coefficients of radial deviation from a perfect ellipse:


where . The result is a set of higher-order terms (, , etc.) describing how the actual isophote deviates from a perfect ellipse. The most interesting coefficients from our point of view are and . The most well-known of these is , which is when the isophotes are “disky” (lemon-shaped) and when the isophotes have a “boxy” shape. The term is non-zero when the boxy/disky shape is rotated with respect to the fitted ellipse; Erwin & Debattista (2013) (hereafter ED13) discuss how both components are affected by the presence of B/P structures in bars.3 We used logarithmic spacing of the semi-major axes for fitted ellipses, to combine high spatial resolution in the central regions and higher at large radii. The ellipse profiles obtained for M31 are shown later in Figures 3 and 4.

3 Method

3.1 Simulations

We want to explore a scenario where the bulge of M31 is a pure B/P bulge or a combination of a classical bulge and a B/P bulge, and in the last case, to constrain the properties of the classical component. That there are no complete analytical descriptions of B/P bulges that reproduce the vertical complexity of these structures, together with our interest in models which evolve in time and naturally develop the B/P bulge structure in equilibrium with a classical bulge, compels us to proceed with N-body simulations. In N-body models B/P bulges emerge from a disk that forms a bar which later buckles vertically, creating a peanut or boxy structure. This is a non-linear process which involves a redistribution of mass in the inner part of the initial disk, where the potentials of all components are involved i.e. the disk, the initial bulge and the dark matter halo. It is not then possible to predict quantitatively the properties of a model after it evolves from its initial conditions. We therefore proceed to make a systematic exploration of the initial parameters with simulations where we change one parameter at the time. We separate here the discussion into the generation of the initial -body models, and the generation of bars and B/P bulges in those models.

-body models: initial conditions

We use the software Nemo, an array of several independent programs and tools to perform -body experiments and analysis in stellar dynamics (Teuben, 1995). To generate the particle models we select the program MaGalie (Boily et al., 2001) based on the method proposed by Hernquist (1993) which solves the Jeans equations to generate galaxies close to dynamical equilibrium with several components, e.g. a bulge, a disk and a dark matter halo. The code works in natural units, ergo the gravitational constant is set to , and the internal units for the variables are: (mass), (distance), (time) and (velocity). The exact values of the scaling factors for converting internal units into physical units vary between each simulation, as explained later in Section 3.2. Typical values are , , and which is . This allows the scaling of the models to M31 by matching velocity and spatial scales independently.

MaGalie builds -body dark matter haloes (DMHs) with different mass density profiles, a cored isothermal profile (used by AB06), or a Hernquist profile. Here we chose the latter, as a convenient approximation to a Navarro-Frenk-White (NFW) DMH profile. It has convergent mass at large radii , and mimics the cuspy NFW profile in the inner parts of the halo as both density profiles behave as within their respective scale radius (Springel et al., 2005). It is given by


where is the radius, is the mass of the halo at and is the scale length. We truncate the haloes at , which defines the actual halo mass in the simulation . The density of the DMHs show some evolution in their outer parts at 20, but quickly stabilises within 100, well before the bar formation, to its final shape.

The initial disk density profile is given by:


where is the cylindrical radius and is the initial radial scale length of the disk, which is fixed by MaGalie to be . The scale height of the disk is and the mass of the disk is also fixed to . As we are interested in the bulge we truncate the disk at . The disks have an exponential radial dispersion profile. We choose a Toomre (Toomre, 1964) value of measured at which avoids axisymmetric instabilities, but allows non-axisymmetric instabilities to grow. We also modified MaGalie to generate and test disks with an initial constant , as explained in Appendix A.

The initial bulges (ICB) are created also with a Hernquist density profile for which, as shown by Hernquist (1990), the projected surface-density profile agrees with a de Vaucouleurs profile (which is a Sérsic profile with index ), within per cent for radii in the range . If integrated, this encloses per cent of the total light. The density profile and parameters are defined here as:


where is the bulge scale length and is the mass of the bulge at . We stop the particle sampling at , which defines the actual ICB mass in the simulation . During the evolution, the ICB density profile near the outer boundary evolves slightly, involving less than 4 per cent of the ICB particles.

-body models: bars & Box/Peanut bulges

Programs like MaGalie can set up models of disk galaxies, but to study the possible coexistence of a B/P bulge with a classical bulge, we need to evolve the initial models to generate the required structures (Athanassoula, 2005). These N-body models generally form a bar that later goes through the buckling instability generating the B/P bulge or thick bar in the centre, which transitions to the thin bar further out that is aligned with the B/P bulge. Transient material trailing or leading the thin bar, like spiral arms attached to the thin bar ends, are not counted as part of the bar. We reserve the term bar for the whole structure, that includes both the thin bar and the B/P bulge.

We want to generate and explore models with bars that show a wide range of boxy structure, pattern speed, bar length and bar strength among others. Therefore, similarly to Bureau & Athanassoula (2005), we choose different concentrations and masses for the DMHs, leading to models dominated by the mass of the disk (MD models), and models dominated by the halo (MH models). MH models usually develop long thin bars and their B/P bulges have a strong X-shape. MD bars are shorter, the thin bar can be very weak, and the B/P bulge has a more boxy shape. We generate and explore B/P bulge models that also include initial classical bulges, as explained in the next section.

-body models: parameter space exploration

We build two sets of models to make a systematic exploration of parameters with a total of 84 simulations. The first set (Set I) contains ICBs combined with B/P bulges and is built from initial models with bulge, disk and DMH components. Here we want to explore how the different ICBs affect the observational parameters. Therefore in this set we vary only the initial mass and size of the ICB component, choosing 12 different masses ranging from 0.05 to 0.6 with steps at every 0.05. For each chosen mass we also explore different sizes for the bulge using 6 values of ranging from 0.1 to 0.35 with steps at every 0.05, ending with a total of 72 simulations for this set. The DMH used in this set has a scale and mass of and .

The second set of models (Set II) contains pure B/P bulges and is built from just disk and DMH initial components. Here we try to generate buckled bars with different boxy structures by changing the concentration and the mass of the DMHs. Therefore we use 3 scale lengths of 10, 15 and 20 and for each we explore 4 different masses : 6, 7, 8 and 9, making a total of 12 simulations.

In addition to these 84 simulations, we have run 100 simulations with different sets of initial parameters, such as disks with different , others with cored isothermal DMHs, and others with an initial disk with initial constant , but found that our fiducial choices best reproduced the bulge of M31, and for conciseness we do not give further details of these simulations here.

Due to the difficulty of plotting the results of the analysis of 84 simulations, we proceed to show only three examples in the next sections, Model 1, Model 2 and Model 3, which belong to Set I and therefore they have the same initial DMH and initial disk. Models 1 to 3 have the same ICB scale length of and differ only in the mass of the ICBs, which are 0.25 (Model 1), 0.05 (Model 2) and 0.5 (Model 3). We also show Model 0 which is a pure B/P bulge of Set II that has the same initial DMH and disk as the previous models. We will show later that Model 1 is our best model for M31’s bulge of all the explored models.

-body models: time integration

To evolve the initial models we used a program also contained in Nemo called gyrfalcON (Dehnen, 2000). Although this program is not parallelized, it is a fast, momentum-conserving tree-code. It uses the same internal units of MaGalie. We choose a time step of and we evolve the initial models until 600 ( 4.65), analysing all the models this standard time. We also analyse and compare some models at 500, 700, 800 and 1000. We choose a tolerance parameter of . For simplicity we use a constant softening parameter of .

The number of particles that we use in both sets for the disk, classical bulge (if present) and DMH are , and respectively, and therefore the respective particle masses for each component are different, with values for Model 1 of , and . To examine the effects of force resolution on our main results we have re-run simulations with new softening parameters using a 50 per cent smaller global and later a 50 per cent larger global . To test the effects of unequal particle masses we follow the prescription by McMillan & Dehnen (2007): the softening for each particle depends on its mass and on the condition of the maximum force () allowed between the particles, obtaining for Model 1 the softenings for the bulge, disk and halo , and . In the resulting simulation with lower resolution we observed no significant variations, while in the simulations with higher resolution we observed that the bar formation is delayed by roughly 100, but the bar evolution, including the buckling, does not change significantly and therefore the age of the bar remains the same, and the results stay unchanged.

3.2 Technique to obtain the best-matching model

In order to compare simulations with observations we also use ellipse on images generated from our simulations. From this we obtain semi-major axis profiles of 5 parameters: the position angle () of the fitted ellipse, the Fourier coefficients and which measure the asymmetry and the boxiness, the ellipticity and the azimuthally average (AZAV) surface-density profile for the simulation images. The position angle is measured anticlockwise with respect to the north celestial pole axis, where . We rely on the coefficient to quantify the strength of the boxiness (negative ) or diskiness (positive ) of the isophotes. As already mentioned in the Introduction, the Sérsic index is useful to classify bulges as classical bulges () or pseudobulges (), and therefore we measured this parameter in M31 and in our models. For this we fit a combination of a Sérsic profile (Sersic, 1968) and an exponential profile to , obtained from M31, and also to from the simulations:


where (Capaccioli, 1989), is the disk scale length, the Sérsic Index, and the effective radius which corresponds to the half light (or half mass) radius of the Sérsic profile. and correspond to the half light and half mass of the Sérsic profile. Here we denote by and the component in Eq.6 fitted by the Sérsic profile. The fit of the parameters is performed with a non-linear least squares (NLLS) minimization method using a Levenberg-Marquardt algorithm, where we explore a full suite of Monte Carlo NLLS realizations with a wide range of initial guesses over all fitted parameters, from which we estimate errors from the standard deviations of the solutions around the best values.

We convert the of the models to surface-brightness dividing by a stellar mass-to-light ratio (). This is determined after the profiles of Eq.6 are fitted to M31 and the models, by scaling of the models to the intensity of M31 measured at the effective radius of M31 () which is , i.e. .

In order to find a best-matching model for M31 we define 6 observational parameters: (1) , corresponding to the difference between the maximum PA () in the boxy region of the bulge and the PA of the disk ; (2) that corresponds to the radius where and the isophotes stop being boxy and start being disky; (3) that quantifies the maximum boxiness of the boxy bulge; (4) the ellipticity at ; (5) the Sérsic index and (6) the effective radius . And finally, we consider an additional parameter, which is the velocity scaling calculated from the line-of-sight maximum velocity dispersion (7) measured in M31.

The procedure used to compare the observations with the models is the following:

  1. We first project the models on the sky as M31, as shown in Fig.1. For this we incline the disk to (where an edge-on disk is , and a face-on disk is ). Then we rotate the projected model around the observer’s line-of-sight axis until the projected disk major axis is aligned with the disk major axis of M31, leaving the position angle of the disk major axis like M31 at (de Vaucouleurs, 1958), and the near side of the disk in the upper part, pointing north-west like M31. We specify the orientation of the model’s bar by an angle in the plane of the disk, such that for the bar is side-on and its major axis is aligned with the projected disk major axis, and increases from the side of the disk major axis at in the direction away from the observer until for the bar major axis is almost aligned with the line-of-sight and is seen nearly end on. Then we generate an image for each model, with a pixel size that slightly varies depending on the model, but with typical values of 5.

  2. We analyse the image of each model with ellipse and measure in the boxy region. is estimated as the error weighted mean of 5 measurements around the maximal PA value (with errors weights from the ellipse fitting), while for the error we use the error weighted standard deviation. This error is larger than the errors estimated by the ellipse fitting, and it takes into account the noise that we observe in the ellipse fitting profiles.

  3. We repeat the previous step for each model, but with different , ranging from 0\degreeuntil 74\degree, until the of each model matches the observed value for M31 which is , obtaining a best bar angle for each model. This parameter is independent of the size scaling which only determines at what distance is located. To determine the error of the best bar angle we calculate where matches the upper and the lower errors of , from which we estimate the error as . As we show later in Fig.3, this error is larger than the effects of the noise in the profile on the bar angle estimation.

  4. We use for each model and we obtain the size scale of each model by matching each to the value for M31. We show later in Section4.2 that the profile of the coefficient can successfully quantify the region where the boxy isophotes of M31 end (), and that our models exhibit a similar behaviour, which makes this value ideal to our interest of restricting the size of the boxy region of M31’s bulge and in our models.

  5. We measure , , and in each model.

  6. We discard the models that do not match the selected 6 observational parameters of M31 (, , , , and ), until we obtain a best model which simultaneously matches the parameters.

  7. We obtain the velocity scale of the best model by matching the maximum value of the line-of-sight dispersion profile along the major axis of the of the model with the value measured from the bulge of M31 by Saglia et al. (2010), i.e. .

  8. We calculate dispersion and velocity profiles and maps for the best model and compare them with M31 observations.

  9. We use the spatial and velocity scaling to obtain the mass scaling and calculate the mass profiles for the best model.

At the end we obtain a model that matches the maximum position angle and the twist of the isophotes, with a boxy region of similar extension and magnitude. And which contains a bulge with similar ellipticity, Sérsic Index, and effective radius. We discuss later the kinematic properties of the best model, which matches the central dispersion and the rotation observed in M31’s bulge.

Figure 1: Schematic diagram of the orientation of the models. We project the models on the sky as M31, giving an inclination to the disk of , locating the near side of the disk pointing to the north-west, and locating also the position angle of the projected major axis of the disk at anticlockwise from the north axis. The bar angle is measured in the plane of the disk. The straight arrow shows the major axis of the bar that is aligned with the projected disk major axis when . The angle increases anticlockwise, as shown by the curved arrow, until for the bar is seen nearly end-on.

4 Results

The results are organized in six subsections. Section 4.1 explains how we match the isophotal twist of the bulge using the parameter , to determine the orientation of the bar in three-dimensional space (). In Section 4.2 we show the results of our morphological analysis with ellipse on the observations and simulations, determining 3 parameters: , , and . In Section 4.3 we use the surface brightness profile of the observations and the surface mass profile to determine two parameters: the Sérsic index and effective radius. Section 4.4 compares the previous mentioned six parameters of the 72 models of Set I, converging in a best matching model. Section 4.5 shows in more detail the properties of the best matching model, showing its kinematics, mass profile and others. Finally in Section 4.6 we analyse the thin bar and the the spurs in the models and compare them with the observations.

4.1 Morphology: bulge isophotal twist & bar angle

Figure 2: Top panel: IRAC 3.6 image and isophotes of M31’s central region. The position angle of the horizontal axis is . The dotted line marks the position of the major axis of the disk (). The blue solid line marks the maximum position angle determined by the fit of the ellipses, reaching , and therefore a difference with the disk of , quantifying the clear isophotal twist in the central boxy region of the bulge. The near side of the disk is located in the upper part of the panel (Walterbos & Kennicutt, 1988). Bottom panel: Image of the central region of the isophotes of Model 1 at 4.65 (600). The mass is converted to luminosity dividing by a stellar mass-to-light ratio of . We reproduce the central twist in the boxy region of the bulge, using the best bar angle of . The blue line marks the maximum of the model which matches . The dashed line marks the projected bar major axis for the best bar angle . The near side of the disk is also located in the upper part of the panel.
Figure 3: Top panel: The projected semi major-axis profiles determined with ellipse for M31 (black), Model 2 (600) (green) with and Model 3 (600) (red) with , with errors in shaded regions. We show two profiles for Model 1 (600), one where (blue dashed-line), and another where (solid blue). The horizontal upper and lower dashed black lines marks the peak of M31 , and the of the disk of . Bottom panel: versus the bar angle . When the bar is seen side-on, and when it is end-on. The horizontal dashed line marks . The dashed curve shows of Eq.7, where the grey area shows the value when . Model 1 values are shown for different (blue up-side triangles). The blue curve is a polynomial fit to estimate the best that matches for Model 1. The blue vertical dash-dot line marks Model 1 best bar angle , where the grey area are the values estimated from the observational errors. The down-side green triangles correspond to Model 2. Dots and crosses correspond to the 72 models of Set I. Dots mark the required for each model to match . The red cross shows the mean and standard deviation for these best bar angle values: . Black crosses are models that cannot reach , due to their concentrated ICBs that generate isophotes with low .

The goal of this step is to match the twist of the isophotes in the boxy region of M31’s bulge relative to the disk, as shown in Fig.2 (top panel). In external galaxies and in -body simulations the presence of a classical bulge, and the orientation of the thin and the thick (or B/P bulge) components of the bar determine the shape and twist of the isophotes of the central region (Bettoni & Galletta 1994, AB06, ED13). ED13 compare surface-density contours of -body barred disk models with isophotes of observed barred disk galaxies using the ellipse fitting analysis and identify several substructures that we also find in our models.

We quantify the twist in M31 and in the models using the difference between the maximum position angle in the boxy region and the disk , obtaining for M31 . This is close to the measurements of Beaton et al. (2007) using the ellipse analysis on the bands J, H and of 2MASS 6X data, obtaining .

We measured the of the iso-density contours in each of our models using different angles for the bar . In Fig.3 (top panel) we show the profile and the for M31 in the boxy region compared to two profiles of Model 1 (0.25) that differ due to the different angles used for the bar, i.e. (blue dash curve) and (blue solid curve). As shown there, this model requires an angle for the bar of to match the observed in M31. We also show Model 3 (0.5) with the more massive ICB which needs an angle of in order to match .

We show in Fig.3 (lower panel) how the bar angle of a 1D-bar measured in the plane of the disk changes its projection into a plane in the sky, with an inclination angle for the disk of , described by the equations:


where is the projection of the bar angle , and is the true position angle of the projected major axis of the bar, that includes the thin bar and the B/P bulge. Therefore if we approximate M31’s bulge as a 1D-bar structure, the required angle to match the photometric twist would be . -body bars and real galaxies are vertically extended and therefore, excluding extreme cases, the difference between the bulge maximum position angle and the position angle of the disk () will usually reach lower values than the infinitesimally thin case (), exactly as we show with our simulations in Fig. 3 (lower panel). In the figure we plot versus for the Models 1 and 2. Using a polynomial fit of order 7 to Model 1 we find that the angle for the bar for which the matches the value observed in M31 is , where the errors are given from the fit using the observational errors, as explained in Section 3.2. This angle generates a twist of the isophotes in the boxy region of Model 1, as shown in Fig.2 (bottom panel). The projected angle of a 1D-bar given by this best angle is , with its error calculated from the average of the upper and lower error in . This locates the of the thin bar (and the B/P bulge) at . Applying the same procedure to Model 2 and Model 3 we recover the best angles for Model 2 and for Model 3. Model 1 and 2 recover more similar values for . Model 3 needs a larger to match , because this model has a more massive ICB that dominates the morphology in the central region, and therefore the isophotes have lower and a less boxy shape and shows a low . Nonetheless, further out it can reach the observed .

Looking carefully at the profiles of the models in Fig.3 (top panel) clarifies why is chosen to determine the bar angle . From all the values in a profile, is the closest value to the estimation given by Eq. 7 (). The deviation of from is shown in Fig.3 (bottom panel). This behaviour is observed in all our models with bars and is the reason why we choose as an indicator for the bar angle. We successfully match the of M31, although the exact radius of the model’s depends on the morphology and length of the thick and thin bar. As a consequence of this choice, the isophotes of the model show a photometric twist slightly weaker than in M31 within the radius where matches . Later, in Section 4.4.2, we show that our conclusions do not change when we increase to produce a more pronounced isophotal twist in the inner part of the bulge region of Model 1.

We repeat this process for all the models, obtaining their respective . The mean and standard deviation of the 72 models of Set I is which shows that the angle does not change much from model to model. Furthermore, we see in Fig.3 (lower panel) that of Eq.7 is a good predictor as a lower limit for bar angles, because none of the 72 models reach values lower than when they match . There are some outliers which never match , reaching always lower values due to the fact that their ICBs have too much mass and/or are too concentrated and their round isophotes dominate.

AB06 used four N-body models with different and compare the spurs generated by the projection of the thin bar of the models with the spur like features at along the major axis of the disk in M31, and concluded that the angle for the bar is between and depending on which model they used. Here instead we use the isophotal twist of the bulge, obtaining , and we argue later in Section 4.6 that structures at are not simply related to the spurs generated by the thin bar.

4.2 Morphology: position angle, boxiness, ellipticity & asymmetry

Figure 4: Parameters determined with ellipse for M31 (black), Model 1 (blue), Model 2 (green), and Model 3 (red) including their errors (shaded areas), within the different regions (vertical lines). The profiles of the models correspond to a snapshot at 600. Top panel: profiles (same as Fig.3) The upper and lower horizontal dashed black lines mark , and . Second panel: profiles. M31 reaches the maximum boxiness at and its boxy region ends at . The models are scaled in order to have the same , but their isophotes may have stronger or weaker boxiness . Third panel: profiles. The upper and lower horizontal lines indicate for the disk and the bulge respectively by Co11 (see Section 4.3) where and . Bottom panel: profiles.

In this section we compare the morphology of the isophotes of M31 with the N-body models. For this we proceed in a similar way to ED13, with a morphological analysis using ellipse on the IRAC 3.6 image shown in Fig.2 and on images of our simulations, obtaining four parameters as a function of the major axis of the fitted ellipses: , , and . The resulting profiles for M31 and the models are shown in Fig.4. We use these parameters to identify isophotal structures and to define five regions (R1 … R5) in M31, going from the inner part of the bulge to the disk. Simultaneously, we plot in the figure the morphological profiles of M31 with three models, Model 1, Model 2 and Model 3 measured at 600, to compare their properties within these regions. These regions and the profiles are similar to the results obtained for M31’s bulge in the bands J, H and K by Beaton et al. (2007), and the morphological properties are also similar to what was found by ED13 for their comparisons between B/P bulges of N-body simulations and other galaxies, using ellipse.

  • Region R1: defined within where and . M31 shows in this region isophotes with a low ellipticity , symmetric, nearly-round isophotes, which makes the uncertain. Model 1 (0.25) shows similar values for , and although the last two parameters indicate slightly more symmetric isophotes than in M31. Model 2 (0.05) has a higher and a more disky shape () in the inner part (), but the boxy region starts sooner (), and is also symmetric . This highlights the necessity of a more massive ICB component, like the one in Model 1, in order to obtain lower and a less boxy structure. But if the mass of the ICB is too high like in Model 3 (0.5) we observe a very low due to the massive ICB, with large fluctuations of the .

    Depending on the scaling parameter of a simulation, region R1 and part of R2 could be within three softening lengths of the centre. Our higher force resolution tests of Model 1 show only small variations of the surface-brightness profile in these regions (see Fig.5) and we find that the isophotes and their morphological parameters change only slightly. The small variation is in part because these are line-of-sight projected quantities and therefore the differences are naturally smaller.

  • Region R2: is the inner boxy region of M31’s bulge, , defined by showing boxy isophotes with , but still roughly symmetric with a low . Model 1 is less boxy than M31 in this region, but equally symmetric and with similar . Model 2 is already more boxy than M31 and its is much larger than observed, reaching already . The models differ again due to the presence of an ICB that dominates in this region, as shown later with the surface density and the mass profiles of Model 1 in Sections 4.3 and 4.5.4. The profiles are better defined, reaching a value of for M31, for Model 1 and for Model 2, but still with some noise. Model 3 still shows a very low and has more disky isophotes compared to M31, Model 1 and Model 2 .

  • Region R3: is the outer boxy region of M31 with , defined by the radius where the isophotes are boxy and start showing an asymmetry () and the radius at where changes from boxy (negative) to disky (positive), indicating the end of the boxy isophotes and the transition to the inner region of the disk, we call this parameter . We also measure in our simulations in internal units and then we use the value of M31 to determine the scaling factor in our models. Therefore all the models end their boxy region at the same place, although they do not necessarily do this in internal units. The maximum boxiness in M31 is located between and .

    Contrary to Region 2 this region shows an asymmetry of the isophotes in M31 given by a positive of , indicating that the is increasing. In our models we see the same features and behaviour of in this region, where two structures overlap which are the outer boxy bulge and the projected thin bar (that generates spurs). Further out the thin bar dilutes into the inner disk, sometimes with transient trailing or leading spiral structures. An increasing (decreasing) indicates that the isophotes are asymmetric and deviating anticlockwise (clockwise) from the major axis of fitted ellipse.

    In this region we can see the twist of the isophotes in M31’s bulge with respect to the isophotes of the disk, as shown in Fig.2. This is reflected in the profile of the ellipses, which increase until a maximum of . In the models we also see that the reaches a maximum in this region.

    The ellipticity of M31 at the effective radius () is , which is in agreement with Co11, where . The profile of Model 1 and agree quite well with the observations in this region. In Section 4.3 we show for this model that the mass of the B/P bulge dominates over the ICB in this region, which has a strong impact on the shape of the isophotes. The profile of Model 2 shows high , where of this model almost doubles the value of M31. Model 3 has generally lower and more disky isophotal shape, and only at starts showing a boxy shape. At its ellipticity reaches the value of M31, but remains almost constant until the disk region.

  • Region R4: This region shows the transition to the disk of M31 at . The isophotes decrease their after reaching . The shows a bump between and which reveals a structure also visible in the isophotes in Fig.2. The models shown here do not reproduce this feature, but we offer some possible explanations later in Section 4.6. The ellipticity of M31 keeps rising, although near the bump it remains roughly constant, and again this is due to the same structure. The ellipticity of Model 2 is higher than M31 and Model 1, reaching already the value of the outer disk . Model 3 has much lower ellipticity than the observations, until .

  • Region R5: the outer part of the disk (), where the reaches , aligning with the line of nodes of the disk. The ellipticity reaches a maximum of (also consistent with Co11, with ). An infinitesimally thin disk would have an ellipticity of using . Our models have a vertically thick disk and reach , as in the observations.

4.3 Photometry: M31’s surface-brightness – two bulge components?

Figure 5: Top panel: The surface-brightness (SB) profiles of M31 (white triangles), Model 1 (600) (blue squares) and Model 0 (600) (cyan circles), within the different regions (vertical lines). We include the inner profile of a re-run of Model 1 with higher force resolution (red circles). The solid curves corresponds to the fit of a Sérsic and an exponential profile, and the dashed curve to the Sérsic component alone. The colours of the curves indicate if the fits belong to M31 (black), Model 1 (blue) or Model 0 (cyan). A stellar mass-to-light ratio of is used in the models to convert into SB. Bottom panel: SB for M31 (white triangles) and Model 1 (blue squares), and each component of the bulge of Model 1 plotted separately: the disk + B/P bulge (magenta solid curve) and the ICB (orange solid curve). We included also the bulge components of Model 1 at the initial time (dashed curves).

In Fig.5 we show the surface-brightness for M31’s bulge region, for Model 0 that is a pure B/P bulge (no ICB), and for Model 1 where we also show its components plotted separately (bottom panel). We also include in the figure our test of Model 1 with a higher force resolution, finding a SB that only differs from the low resolution model by 0.5 per cent within 20. We determined the Sérsic index and the effective radius for the models and for M31. We convert the surface density of Model 1 to SB dividing by a stellar mass-to-light ratio calculated as in Section 3.2, (), which agrees with the range of values estimated by Meidt et al. (2014) for the IRAC 3.6 band, which goes from to for a Chabrier IMF, depending on the metallicity and the age of the stellar population. The mass-to-light ratio shifts the profile of the models in the vertical axis of Fig.5 and is an important value, but not critical at this step of the analysis, as it does not change the shape of the profile.

We see in the figure that the SB of Model 1 agrees quite well with M31 within Region 3 (), where the B/P bulge dominates. Model 0 also shows a similar SB. Within Region 2 () M31 keeps increasing its SB, while the pure B/P bulge, Model 0, does not. Model 1 matches the observed SB due to the contribution of its ICB. Within Region 1 () the SB of M31 keeps rising and it is matched by Model 1, where the ICB component now dominates with its cuspy SB over the B/P bulge component that shows a cored SB profile. It’s interesting that Region 1 was defined with a morphological analysis of the isophotes of M31, and yet we see that this region is exactly where the SB of the classical bulge component dominates in Model 1. Model 0 shows also a cored SB and cannot reach the SB of M31, which highlights the importance of the ICB component and the mass concentration of this component. We also show the SB of the two bulge components of Model 1 at the initial time, before the formation of the bar. The mass redistribution of the initial disk when its material forms the bar increases its SB in the central regions by . The ICB almost does not change its SB despite the dynamical events like the bar formation and buckling, and the subsequent angular momentum transfer between components.

We also show in Fig.5 the profiles that we fit to the SB of M31, Model 1 and Model 0. From the fit to M31 we obtain the Sérsic index , the effective radius , and the disk scale length of . and converted to surface-brightness are and . These values are consistent with the results of Co11, where they used several fitting methods on surface-brightness profiles generated with different fields or cuts of the image, and depending on the masking and the fitting method, finding values, varying from to , from to and a range of disk scale lengths of to .

We find in our models that after the bar formation the surface-brightness profiles in the disk region shows a broken exponential profile, with two disk scale lengths. In the outer region () we find a that is similar to the initial one, while in the inner region () we find a larger . This is also present and explained by Debattista et al. (2006) as part of the secular evolution due to the mass re-distribution of the disk after the bar formation. Therefore, as we seek models of the bulge region of M31, we choose to leave the disk scale length as a fixed quantity, using , and obtain for M31 a Sérsic index and an effective radius . We then fit the models in the same way, leaving the parameters of the Sérsic profile and the central disk intensity () free to vary. In this way we can compare the parameters of the Sérsic component of observations with the simulations alone, which is the component that better quantifies the properties of the bulge’s density profile.

We find that Model 0, without a ICB component, has a Sérsic index , lower than M31. On the other hand Model 1, with its ICB, behaves very similarly to M31, with and . We find that models of Set I that have higher and smaller than Model 1 tend to show higher than observed. On the other hand, models of Set II, which are pure B/P bulges rarely show higher than 1 (Debattista et al., 2006), and therefore this set is ruled out. In the outer parts, beyond Region 3, Model 1 starts to show lower densities than M31’s disk. We find that our N-body disk models tend to decrease faster than the observed profile in M31. This disagreement is important, but not critical as our interest in this paper is to find models for M31’s bulge.

We also find that the Sérsic index of a model can change with the orientation of the bar. When the bar is side-on () the Sérsic index can be higher, because the ICB has less material of the B/P bulge in the line of sight, emphasizing its steep density profile. And when the bar is end-on () the Sérsic index can be lower, because more material of the B/P bulge is in the line of sight, which has roughly an exponential () profile. Furthermore, the Sérsic index decreases also in time as the bar and B/P bulge grows incorporating more material, hiding the ICB component (depending also on the orientation).

4.4 Parameter space for the ICB

Selection of the best model

Figure 6: Model - M31 comparison in the parameter space of the ICB versus size , for the 72 models of Set I at 600. All models already match and (except for 6 cases shown with a purple filled circle in the upper left corner). Blue, brown, magenta and red represent the parameters , , and . The parameters of the models that agree well with the parameters of M31 are shown with filled coloured circles, and they have values that are within the range . Solid thin circles mark models with values larger than the parameters in M31 by , while dashed circles mark values lower than M31 by . Model 1 matches simultaneously all parameters and is the best matching model (square).

As explained in Section 4.3, the pure B/P bulges do not show high mass concentrations in the centre and therefore their Sérsic indexes are much lower than observed, all with values , ruling out all models of Set II. Now we investigate the properties of the ICBs of Set I that can simultaneously reproduce all the morphological observables of M31’s bulge. The most important initial parameters of the ICB are its mass and its size . We show these two parameters in Fig.6 for the 72 models of Set I. The different colours used in the figure identify the four parameters: , , and the Sérsic index . In order to compare the models and M31 we use the fractional difference between data and model (in percent): (where corresponds to the each of the four parameters). If the circle is fully coloured with one colour it means that this parameter agrees to within of the value in M31. Empty circles marked by a solid line (or dashed line) denotes that this parameter is larger than in M31 by (or lower than in M31 by ). The threshold of 15 per cent is chosen for two reasons: (i) typical errors in the parameters are between 10 and 20 per cent, and (ii) it is roughly at this value where is found a unique candidate that simultaneously matches all the parameters. By construction the models of Fig.6 already match and (except for a few models marked in the figure with purple circles that have ICBs that are too concentrated to match , as explained Section 4.1). We highlight with a square in the figure the best model, Model 1, which simultaneously reproduces all the parameters.

The Sérsic index gives information about the mass concentration of the models. We explore a reasonable range of sizes and masses for the ICB, which give as low as 1.13 and as high as 3.23, with an average and standard deviation of . As expected the ICB with sizes smaller than are too concentrated, resulting in too high, with values between 2.52 and 3.23. Fisher & Drory (2008) showed that pseudobulges and classical bulges can be distinguished by the threshold . The pure B/P bulge models fail to reproduce M31 surface-brightness, indicating that a classical bulge component is needed, but its mass contribution is limited in order to match .

The boxiness tends to be weaker in the upper right corner of the figure (), where the ICB are more massive and concentrated, which makes the isophotes less boxy. Also, a compact bulge can help to stabilize the buckling instability, therefore making B/P bulge less boxy (Sotnikova & Rodionov, 2005). The inverse is also true, as shown by the models with less concentrated ICB and with more strongly boxy structure (or even X-shape) that are located in the bottom left corner of the diagram (). The ellipticity behaves similarly to , where more compact ICB give round isophotes and therefore low . We find that the variations of the parameters , and vary more smoothly in the diagram, contrary to which is more scattered. The boxy structure forms in a non-linear process and its size can show a wide range of values of in internal units (), which in the case of Set I varies from 1.26 to 3.74. Since is rescaled along with , this results in a wide range for as well.

Best model parameters & bar angle

In the previous section we analysed and compared the 72 models of Set I with parameters of M31, and at the end of the selection process only one model remains, Model 1. We found that the bar angle