Andromeda chained to the Box  Dynamical Models for M31: Bulge & Bar
Abstract
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 Nbody 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 threedimensional halfmass radius . The B/P bulge component has a mass of and a halfmass 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 10ring 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 .
keywords:
galaxies: spiral — galaxies: bulges — galaxies: bar — galaxies: kinematics and dynamics — galaxies: individual (Andromeda, M31, NGC224)) — galaxies: Local Group – galaxies: structure — methods: nbody1 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 nonrotating 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 Nbody 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 surfacebrightness 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 bulgetototal 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 Nbody 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 setup 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 largescale 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 northwest, 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 semimajor axis , this is done by iteratively adjusting the ellipse centre (pixel position in and ), orientation (position angle ), and ellipticity ( or , where is the semiminor 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 surfacebrightness profile (SB) using the 3.6 zeropoint 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:
(1) 
where is the mean intensity along the bestfit ellipse; the and components are automatically zero for the bestfit ellipse. In fact, the actual computation divides the and higher coefficients by the local radial intensity gradient and the ellipse semimajor axis to generate coefficients of radial deviation from a perfect ellipse:
(2) 
where . The result is a set of higherorder 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 wellknown of these is
, which is when the isophotes are “disky” (lemonshaped)
and when the isophotes have a “boxy” shape. The term is
nonzero 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 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 Nbody simulations. In Nbody 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 nonlinear 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 NavarroFrenkWhite (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
(3) 
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:
(4) 
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 nonaxisymmetric 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 surfacedensity 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:
(5) 
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 Nbody 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 Xshape. 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, momentumconserving treecode. 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 rerun 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 bestmatching model
In order to compare simulations with observations we also use ellipse on images generated from our simulations. From this we obtain semimajor 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) surfacedensity 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:
(6)  
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 nonlinear least squares (NLLS) minimization method using a LevenbergMarquardt 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 surfacebrightness dividing by a stellar masstolight 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 bestmatching 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 lineofsight maximum velocity dispersion (7) measured in M31.
The procedure used to compare the observations with the models is the following:

We first project the models on the sky as M31, as shown in Fig.1. For this we incline the disk to (where an edgeon disk is , and a faceon disk is ). Then we rotate the projected model around the observer’s lineofsight 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 northwest 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 sideon 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 lineofsight 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.

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.

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.

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.

We measure , , and in each model.

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.

We obtain the velocity scale of the best model by matching the maximum value of the lineofsight 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. .

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

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.
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 threedimensional 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
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 surfacedensity 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 isodensity 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 1Dbar 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:
(7)  
(8) 
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 1Dbar 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 1Dbar 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 Nbody 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
In this section we compare the morphology of the isophotes of M31 with the Nbody 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 Nbody simulations and other galaxies, using ellipse.

Region R1: defined within where and . M31 shows in this region isophotes with a low ellipticity , symmetric, nearlyround 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 surfacebrightness 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 lineofsight 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 surfacebrightness – two bulge components?
In Fig.5 we show the surfacebrightness 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 masstolight 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 masstolight 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 surfacebrightness are and . These values are consistent with the results of Co11, where they used several fitting methods on surfacebrightness 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 surfacebrightness 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 redistribution 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 Nbody 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 sideon () 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 endon () 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
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 surfacebrightness, 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 Xshape) 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 nonlinear 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