Unravelling stellar populations in the Andromeda Galaxy

To understand the history and formation mechanisms of galaxies it is crucial to determine their current multidimensional structure. Here we focus on stellar population properties, such as metallicity and [$\alpha$/Fe] enhancement. We devise a new technique to recover the distribution of these parameters using spatially resolved, line-of-sight averaged data. Our chemodynamical method is based on the made-to-measure (M2M) framework and results in an $N$-body model for the abundance distribution. We test our method on a mock dataset and find that the radial and azimuthal profiles are well-recovered, however only the overall shape of the vertical profile matches the true profile. We apply our procedure to spatially resolved maps of mean [Z/H] and [$\alpha$/Fe] for the Andromeda Galaxy, using an earlier barred dynamical model of M31. We found that the metallicity is enhanced along the bar, with possible maxima at the ansae. In the edge-on view the [Z/H] distribution has an X shape due to the boxy/peanut bulge; the average vertical metallicity gradient is equal to $-0.132\pm0.006$ dex/kpc. We identify a metallicity-enhanced ring around the bar which also has relatively lower [$\alpha$/Fe]. The highest [$\alpha$/Fe] is found in the centre, due to the classical bulge. Away from the centre, the $\alpha$-overabundance in the bar region increases with height, which could be an indication of a thick disc. We argue that the galaxy assembly resulted in a sharp peak of metallicity in the central few hundred parsecs and a more gentle negative gradient in the remaining disc, but no [$\alpha$/Fe] gradient. The formation of the bar lead to the re-arrangement of the [Z/H] distribution, causing a flat gradient along the bar. Subsequent star formation close to the bar ends may have produced the metallicity enhancements at the ansae and the [Z/H] enhanced lower-$\alpha$ ring.


Introduction
The central parts of disc galaxies are occupied by bulges, which can be classified into two broad categories (Kormendy & Kennicutt 2004;Fisher & Drory 2016). The classical bulges were probably formed very early on, from violent early gas-rich mergers or mergers within clumpy disks (Hopkins et al. 2009;Brooks & Christensen 2016;Bournaud 2016). Boxy-peanut and disky bulges are thought to be built through evolution of the disc component, triggered by bar formation (Kormendy 2013;Fragkoudi et al. 2020). It has been found that different types of bulges can coexist in a single galaxy (e.g. Erwin et al. 2015). In such a case, the bar will transfer some of its angular momentum and spin-up the classical bulge (Saha et al. 2012(Saha et al. , 2016. A major channel of bar formation is a global instability. N-body simulations showed early-on that disc galaxies are prone to development of elongated structures in their centers (Miller et al. 1970;Hohl 1971). Shortly after their formation bars thicken, acquiring a boxy/peanut (b/p) shape in the sideon view (Combes & Sanders 1981;Combes et al. 1990). This is a result of another instability called buckling (Raha et al. 1991;Athanassoula & Misiriotis 2002;Debattista et al. 2006). An alternative explanation of this process is thickening through resonances (Combes et al. 1990;Quillen et al. 2014). The verti-⋆ e-mail: ggajda@mpe.mpg.de cally extended part of the bar constitutes the b/p bulge (see , for a recent review). More about the theoretical understanding of bar physics can be found in the reviews by Athanassoula (2013) and Sellwood (2014). The fraction of barred galaxies grow with cosmic time, starting from 10% at z = 1 (Sheth et al. 2008;Melvin et al. 2014) to about 50 − 70% in the local Universe (Skibba et al. 2012;Erwin 2018). Also the abundance of the b/p bulges grows with time (Kruk et al. 2019), reaching ∼ 40% (see also Lütticke et al. 2000).
The discussion about the impact of bars on the stellar populations has not yet concluded. Some of the authors suggest that bars lead to higher metallicities in the galaxy centers (Moorthy & Holtzman 2006;, while the others did not identify significant differences (Jablonka et al. 2007;Williams et al. 2012;Cheung et al. 2015). Pérez et al. (2009) found all types of metallicity gradients along the bars: positive, flat and negative. Williams et al. (2012) argued that the gradients in the bars are flatter than in the discs, suggesting that the bar formation smears out pre-existing gradients. Coelho & Gadotti (2011) concluded that the bulges in the barred galaxies are on average 4 Gyr younger than in the unbarred ones. Sánchez-Blázquez et al. (2014) did not find any differences in the metallicity gradients in the outer parts of barred and unbarred galaxies, contrary to some simulation predictions (Friedli et al. 1994;Minchev & Famaey 2010). In the Milky Way, stars in the immediately surrounding disc are a bit more metal rich than in the bar region (e.g. Hayden et al. 2015). A&A proofs: manuscript no. unravelling_m31 Our neighbour, the Andromeda Galaxy, is an excellent target for investigating the stellar populations in a center of a large galaxy. M31 has long been described as hosting a classical bulge (e.g. Kormendy & Bender 1999;Kormendy et al. 2010). However already early-on, Lindblad (1956) posed that the central twist of the M31's isophotes is caused by a bar. This argument was strengthened by Athanassoula & Beaton (2006), who compared barred N-body models to the near-infrared image of Beaton et al. (2007) and concluded that the bar has length of ≈ 1300 ′′ (∼ 5 kpc). Blaña Díaz et al. (2017) considered an array of models consist of both a classical bulge and a bar. They concluded that the classical bulge contributes ∼ 1/3 mass of the Andromeda central part, while the b/p bulge contribution is ∼ 2/3. Blaña Díaz et al. (2018) extended that work, fitting both the Barmby et al. (2006) infrared image and the kinematics derived by Opitsch et al. (2018), using the made-to-measure (M2M) technique. They concluded that the bar has length of ≈ 4 kpc and is oriented at 54.7 ± 3.8 • with respect to the line of nodes of M31. Opitsch et al. (2018) provides a more extensive account of the evidence for the barred nature of M31. Several lines of evidence, such as presence of the Giant Stream (e.g. Hammer et al. 2018), the recent burst of the star formation  and the stellar age-velocity dispersion relation (e.g. Bhattacharya et al. 2019) point to a recent (∼ 3 Gyr ago) merger with a mass ratio approximately 1:5, which could have also left an impact on the distribution of the stellar populations in the center of M31.
Recently, a wide IFS survey of Andromeda was performed by Opitsch et al. (2018). Subsequently, Saglia et al. (2018) analysed their spectra to study stellar populations in M31. They found that 80% of their measurements had ages larger than 10 Gyr. The metallicity along the bar was solar with a peak of 0.35 dex in the very center. The [α/Fe] enhancement was approximately 0.25 dex everywhere, rising to 0.35 dex in the center. They proposed a two-phase formation scenario, according to which at first the classical bulge and proto-disk formed in a quasi-monolithic collapse or a violent instability. Later, the bar formed and buckled into a b/p bulge and the star formation continued not only in the disc, but also in the inner 2 kpc.
Galaxies are distant objects and we can observe them only in projection on the sky. However, to really understand their structure we need to decipher the three-dimensional distribution of their components. It has been demonstrated that a deprojection of the surface density is increasingly degenerate away from special cases such as a thin disk or an exactly edge-on axisymmetric system (Rybicki 1987;Gerhard & Binney 1996). For triaxial systems already Stark (1977) illustrated the degeneracy by finding a sequence of ellipsoidal bulge models that would reproduce the observed twist between the bulge and the disk isophotes in M31, given a common principal plane. Besides the density distribution of the luminous and dark components, the distribution of the stellar population properties is also of interest. In particular, metallicity, elemental abundances and stellar ages are vital to understanding the evolution of galaxies.
Here we consider the determination of the three-dimensional distribution of mean stellar population properties from the observational data for M31. We use the made-to-measure technique (Syer & Tremaine 1996) to incorporate the constraint that in dynamical equilibrium stellar population properties must be constant along orbits. With this method one adjusts a particle model to a set of constraints by iteratively adjusting the masses of the particles. The technique was adapted by de Lorenzi et al. (2007) to fit observational data through minimization of a respective χ 2 and implemented as the nmagic code.
It has been used to study elliptical galaxies (de Lorenzi et al. 2008;Das et al. 2011), Milky Way (Portail et al. 2015, 2017a and Andromeda (Blaña Díaz et al. 2018). The M2M method was found by Long & Mao (2012) to give similar results to the Schwarzschild (1979) modelling and to reproduce well analytically known distribution functions (Tagawa et al. 2016). In particular, Portail et al. (2017b) used distance resolved stellar parameter data to reconstruct the distribution of metallicities in the Milky Way bulge. The Schwarzschild orbit method has recently been extended to stellar population modelling as well (Poci et al. 2019;Zhu et al. 2020).
We build on the M31 dynamical model by Blaña Díaz et al. (2018) and the stellar population analysis by Saglia et al. (2018). Our goal is to construct a three-dimensional model of metallicity and α-enhancement in the Andromeda Galaxy. In Section 2 we present and discuss the available data based on Lick indices. Next, in Section 3, we present our newly developed technique, test it on mock galaxy model and comment on the uncertainties introduced by the limited nature of the available information. Then we apply our method to M31, first to [Fe/H] in Section 4.1 and then to [α/Fe] in Section 4.2. We discuss our results and plausible origins of the observed trends in Section 5, and finally summarize our conclusions in Section 6.

Spatially resolved stellar population maps for M31
The Andromeda Galaxy is the closest big spiral, which is both an opportunity and a challenge. The close distance enables us to create very detailed maps of various quantities. On the other hand, M31 has a large size on the sky, thus one needs multiple visits or a survey with an extended sky coverage. As the constraints for our model we use the publicly available stellar population properties of the Andromeda Galaxy derived by Saglia et al. (2018), based on the data collected by Opitsch et al. (2018), using the W-VIRUS instrument (Fabricius et al. 2012). They covered mostly the bulge area and sparsely sampled the disc along six directions. Spectra were rebinned to reach minimum S/N = 30 and the analysis yielded 6473 usable spaxels. Saglia et al. (2018) measured the absorption line strengths in the Lick/IDS system (Worthey et al. 1994 Trager & Somerville (2009) found that metallicities derived from the Lick indices, so-called SSP-equivalents, follow the stellar mass-weighted metallicity of the given spectrum. Metallicity obtained in this way slightly underestimates the true value and has scatter smaller than 0.1 dex.
While it might not be surprising that one can treat the measured [Z/H] abundance as mass-weighted averages of the underlying stellar populations, certainly the [α/Fe] abundance ratio needs more explanation. Firstly, Serra & Trager (2007) found that [α/Fe] measured from the Lick indices well-reproduces a light-weighted mean of the stellar populations. Secondly, Pipino et al. (2006Pipino et al. ( , 2008 argued that in the case of the αenhancement the light-and mass-weighted averages should give the same results. On the other hand, Trager & Somerville (2009) found that ages derived from the indices can be severely underestimated (compared to the mass-weighted mean) in the presence of young stellar populations. Moreover, 40% of Saglia et al. (2018) age measurements fall on the edge of the grid of models at 15 Gyr. As a result, the age map is mostly featureless, except for the dust influence. Taking everything into consideration, we decided against fitting the age distribution of M31. MacArthur (2005) concluded that the impact of the dust on a single measurement is probably smaller than the typical statistical uncertainty, but here we are going to use the measurements in conjunction. Thus, the uncertainties we are going to report here will be smaller than the ones of MacArthur (2005), hence the dust-related bias may have an impact on our results.
Additionally, the Lick indices suffer from the well-know agemetallicity degeneracy, which may bias the [Z/H] measurements in a way that is not accounted for using simple gaussian uncertainties. The issues with the retrieval of the stellar population parameters could be more successfully addressed using the full spectral fitting (MacArthur et al. 2009;Sánchez-Blázquez et al. 2011;Cid Fernandes et al. 2013, 2014 The Saglia et al. (2018) statistical uncertainties were estimated from the relevant χ 2 distributions, separately upper and lower limits. Initially, as a statistical uncertainty we took the larger of the two. We calculated also a local uncertainty, i.e. for each cell we computed the standard deviation of the distribution of its neighbours. Finally, we derived also an asymmetry uncertainty. For each pixel of each cell located at (R x , R y ) we found what is the value of the parameter at (−R x , −R y ) (if it exists). We averaged those two values and took the half of the difference between the given cell and its reflection as the asymmetry estimate. In the end, as the final uncertainty we took the largest of the three. Saglia et al. (2018) made available the original positions of the W-VIRUS fibres and their allocation to the voronoi cells. Since we needed to divide the sky plane into cells corresponding to the spectra, we came up with the following procedure. First, we divided the plane of the sky into a fine grid of square pixels of 1 ′′ size. Then, for each pixel we found what is the closest observed fibre, provided that it was closer than 5.35 ′′ . This value ensures that all pixels inside the triangular fiber pattern of W-VIRUS are properly assigned. If for a given pixel the closest fibre is farther away, we treated that pixel as missing and we do not further use it in our considerations.

M2M modelling
Our aim in this contribution is to construct an N-body model of a stellar populations parameter φ, for example metallicity, α-enhancement or age. To achieve this, we are going to use an observed map of the mean of the given parameter. Such a measurement is believed to be robustly obtained from full spec-tral fitting (e.g. Cid Fernandes et al. 2013, 2014 or Lick indices (Trager & Somerville 2009). For each line of sight j (e.g. a pixel or a voronoi cell) we are going to use a mass-weighted mean Φ D j and its associated uncertainty σ Φ, j .
Furthermore, let us assume that we already have a dynamical, equilibrium N-body model of the galaxy, obtained through fitting the surface density and the kinematic data. We are keeping the dynamical model fixed and we are not using the stellar populations to alter the model. In other words, in the M2M context, we are basically keeping the particle mass weights constant. This is an important point, since we will require our stellar population model to be consistent with the orbital distribution, therefore we will partially lift some of the degeneracy related to a deprojection of an image into a fully three-dimensional distribution. Now let us assign to every particle i a single value of the parameter of interest φ i (e.g. [Z/H], [α/Fe] or age). We then observe our model galaxy from the same distance, at the same viewing angles and through the same lines of sight as the real galaxy we are considering. For each line of sight j we calculate the model observables as where m i are the masses of the particles and the sums are performed over all the particles present along the line of sight j. Such a measurement might be quite noisy, thus we replace it by a time-averaged value, as originally proposed by Syer & Tremaine (1996) where τ is a constant chosen in relation to the dynamical timescale. In practice, we approximate the integral by a discrete rule, which updates the value after each iteration (see Syer & Tremaine 1996;de Lorenzi et al. 2007). The value "observed" in the model is then compared with the data using a χ 2 statistic where the summation goes over all of the observed lines of sight.
In the usual made-to-measure manner, we want to change the particle values of the parameter φ i so that they fit the data optimally. We construct a merit function F = − 1 2 χ 2 and while the particles orbit in the galactic potential, we apply the following force-of-change where ǫ is a suitably-chosen numerical parameter. In practice, we apply this equation, using the Euler method, in regularly spaced intervals, which we call iterations. In the dynamical formulation of M2M it is common to use an entropy term in the merit function (Syer & Tremaine 1996) in order to reduce the width of the particle mass distribution. However, it is natural to expect a non-negligible width and skewness of the distribution of stellar properties. Hence, we decided not to include any entropy term.
The possible values of φ i that a particle can have should be limited for both physical reasons and due to limitations of the A&A proofs: manuscript no. unravelling_m31 technique used to obtain the data. For example, one would limit ages to values of 0 -13.8 Gyr, while metallicity and abundance would be limited by the extent of the stellar library.
As in Blaña Díaz et al. (2018), we calculate the potential using the hybrid method of Sellwood (2003). We combine a polar-grid solver of Sellwood & Valluri (1997) (updated by Portail et al. 2017a to accommodate different softening lengths in the radial and vertical directions) and a spherical harmonic solver of de Lorenzi et al. (2007). We let the potential rotate around the minor axis of our model with the angular velocity equal to the pattern speed of the bar. In such a potential we evolve the positions and the velocities of the particles, effectively treating them as test particles.
The fitting procedure follows a usual M2M route. First, we initialize all of particle φ i . Then, we let the model evolve for N smooth iterations, so that the observables are properly smoothed. Next, for N fit iterations we fit φ i of the particles, following (4). Finally, we let the model relax for N relax iterations, so we can check if it was not overfitted. From the final values of the particle φ i we compute other interesting characteristics, such as profiles and deprojected maps.
Our method can be summarized as follows. We start with an N-body model that is a faithful representation of a galaxy. We tag every particle with a single value of e.g. metallicity. As the particles move in the galactic potential, we adjust their metallicities to fit the observed map of the mean metallicity in the galaxy.
Initially, we made a number of tests with moderately inclined mock galaxies (i.e. not edge-on). We found that if we used a very simple initialization of φ i , making it equal to a constant value everywhere, our technique was not able to recover on its own the correct vertical gradient of φ. To overcome this issue, we initialize the values of particle φ i as a function of its height above the plane of galaxy where G and N are constants, z i is the particle's vertical coordinate and z 0 is a normalisation constant, equal in our case to the mean absolute vertical coordinate of all particles (which is equal to the scaleheight in the case of the exponential profile). We try a set of possible G and N and check which one results in the lowest value of final χ 2 . Then we use this initial condition for final results and uncertainty estimation. One could wonder if a linear function of the vertical coordinate is sufficient, or should we use a different, possibly more complicated function. Unfortunately, we are not aware of any observationally or theoretically motivated functional form for the vertical profiles of metallicity or α-enhancement. We decided to use the next simplest polynomial (after a constant value) of degree one, which has two free parameters. As we validated in Sec. 3.2, in this way we are able to capture most of the variation, but not small details. Another scheme may be better, but one would need more data to judge it. We note that the prior is not a function of the radial coordinate. We tested what would change if we had an initial condition with some radial gradient. It turned out that the radial profile is very well constrained by the data and such a prior would not matter.
To estimate the uncertainties of e.g. profiles, we use the following procedure. We initialize particle φ i using the best-fit vertical profile with an additional gaussian noise to seed randomness. To the data values Φ D j we add a gaussian noise with zero mean and standard deviation of σ Φ, j . Next, we refit this new data and recompute the values of interest, e.g. profiles. We repeat this procedure 100 times and from the variance of the profiles we estimate its uncertainty. To such a statistical uncertainty we add in quadrature a spread of the profiles that were obtained from models with different initial vertical profiles and were within 1σ from the χ 2 minimum.

Tests on mock data
When a new method is proposed it should be verified on suitable mock data, so one can be reasonably convinced that it gives correct answers. Hence, here we describe our tests.
To create mocks we used the chemodynamical barred galaxy model created by Portail et al. (2017b). It consists of 10 6 stellar particles in dynamical equilibrium with its dark matter halo. Originally, it was a disc galaxy with a bar of 5 kpc length. Since we wanted to make a comparison to M31, which has a 4 kpc bar, we decided to adjust the extent of the model. We scaled all of the sizes by a factor of 4/5, all of the velocities by also 4/5 and masses by (4/5) 3 1 . We observed the model at the distance (785 kpc), the inclination (77 • ) and the position angle of the bar with respect to the line of nodes (54.7 • ) the same as in M31.
In Portail et al. (2017b) each particle has four weights, representing fractions of the particle mass corresponding to four bins of the stellar [Fe/H]. For the purpose of this test we assign each particle a single mean metallicity that reflects the fractional weights. Thereby, we obtain a reasonable model of mean metallicity in a barred galaxy. We use the same set of the voronoi lineof-sight cells as Saglia et al. (2018) and we let the model evolve for 10 4 it (1 it = 1.1 × 10 −4 Gyr) to smooth the observables. In order to create a realistic observed [Fe/H] map we added a gaussian noise with zero mean and σ [Fe/H] = 0.04 dex, equal to the average uncertainty of [Z/H] reported by Saglia et al. (2018).
We limit the possible values of [Fe/H] to the same range as the [Z/H] grid of models in Saglia et al. (2018), i.e. from −2.25 to 0.67 dex. As the underlying dynamical model we use the model from Portail et al. (2017b), thus we do not have any additional uncertainty that would arise if the model used in the fitting did not correspond to the density distribution of the "data".
We illustrate our procedure of fitting the vertical profile in Figure 1. We start with imposing the vertical metallicity pro-1 Recall that in dynamics there are three basic dimensions, which can be chosen as e.g. length, velocity and mass. The gravitational interactions are invariant under the transformation x → αx, v → βv and m → αβ 2 m, where x denotes coordinates, v denotes velocities and m denotes masses. Note that due to this transformation time t → (α/β)t. files, showed by the dashed lines. It may seem that they are far away from the vertical metallicity profile of the original mock galaxy (black solid line). However, in the initial smoothing phase the gradients get shallower, due to the phase mixing of the particles. Moreover, the prior is important mostly where the data constraints are weak, i.e. at large heights and large distances from the centre. After running the code, each initial profile results in a slightly different final one (solid lines) and a different final value of χ 2 . In the Figure 1 we show models with different initial gradient G, but the same normalisation N at z 0 = |z| = 0.297 kpc. The initial linear profiles are transformed into more complicated functions. The best fit (i.e. the lowest χ 2 ) is obtained for an initial G = −0.18 dex/kpc, while the other two correspond to 1σ worse cases. In fact, the blue line approximates what happens if we use an initialization equal to a constant value. While at a high inclination such an initialization is not much worse, we reiterate that it is noticeably worse at more face-on cases. Similarly, we find the best-fit value of N = 0.21 ± 0.03 dex. The largest difference between the best-fit profile and the other models inside 1σ region is considered as a part of the final uncertainty. We want to stress again that the exact values of G and N are not directly related to the final vertical profile. We checked the impact of the vertical prior on the radial profile and it turned out that in the central part it is rather minor. In the outer part it is more noticeable, which is reflected by a wider uncertainty band for the radial profile beyond R ∼ 7 kpc.
In Figure 2 we present both the mock data and the fitted map of mean metallicity. Reduced χ 2 red = 0.99 indicates an almost perfect fit. However, such a good result is not surprising since the dynamical model is correct and uncertainties of the data are perfectly gaussian, with known amplitudes. It seems that our proce- We took into account all stellar particles with |z| < 3 kpc. In the second panel we plot the vertical profile for all particles R < 5 kpc. In the last panel we present the azimuthal profiles for the particles with R < 5 kpc and |z| < 3 kpc.
Our fitting routine is able to recover all of the features of the radial profile, including [Fe/H] peaks at ≈ 1 kpc and close to the end of the bar at ≈ 5 kpc, as well as further negative gradient up to 10 kpc. This could be surprising, since the data along the major axis reach only 5 kpc, however the key lies in the fields along the minor axis, which reach large distances due to a projection A&A proofs: manuscript no. unravelling_m31 effect at the inclination of 77 • . In the azimuthal profile, the main variation is related to the enhancement along the bar, which is aligned with ϕ = 0, ±π.
As we could already see in Figure 1, the vertical profile is not reproduced exactly. However, we are able to recover an average vertical gradient. One could also infer certain degree of non-linearity in the profile, which is flatter closer to the galaxy plane than at larger heights. We made another test to judge if the vertical profile is unduly influenced by a limited on-sky data coverage. We constructed a complete dataset that is similar to Fig. 2, but without any white spots whatsoever, i.e. filling the whole 6 by 12 kpc field. However, for technical reasons we used four times less voronoi cells. The yellow lines in Fig. 3 depict the result of such a fit. The vertical profile profile is closer to the original model, especially at low heights. A discrepancy remains, which could be ascribed to other reasons, such as the orbital degeneracy discussed in the next subsection, the small spatial resolution or the noise. It may be surprising that the recovery of the vertical profile in case of an almost (77 • ) edge-on galaxy is problematic. However, recall that our early tests showed that the constraints on the vertical profiles in more face-on galaxies are much weaker. Secondly, the vertical variation leaves rather small imprint on an on-sky map, which is dwarfed by the radial changes.
The way we calculate the uncertainties also seems to be justified, judging from the radial and azimuthal profiles. On the other hand, the uncertainty in the vertical profile at small heights appears to be underestimated. We want to note that the uncertainties in a single profile and between different profiles are highly correlated, since we are using spatially-extended orbits as building block of our M2M model.

Discussion of the method
One could ask a question what we can actually constrain with this method, using maps of mean Φ D value on the sky. First, let us remind how one could understand M2M dynamical models, which we use as a basis for our work. They are obviously a set of particles in a dynamical equilibrium. However, they can be also regarded as a collection of orbits. Each orbit has an associated mass (equal to the mass of the particle) and the sum of the orbital masses gives a density distribution, whose potential generates the aforementioned orbits. This interpretation resembles the Schwarzschild (1979) modelling, which was used by Long & Mao (2018), Poci et al. (2019) and Zhu et al. (2020) to study stellar populations. We note that Long (2016) implemented a similar technique, however it was applied to absorption line strengths, limited to symmetrized data and not tested besides converging to an acceptable χ 2 .
If we regard our system as a collection of orbits {i}, then we are actually fitting the mean values of that parameter on a given orbit. However, we cannot constrain a distribution function of the parameter on the orbit, e.g. if it is δ-like or very broad. Therefore, we also cannot derive the distribution function of φ for the whole galaxy, we can only give a rough lower limit on its width. As an extreme example, one can imagine that a spatial map of a galaxy is constant everywhere and equal to φ 0 . One would concluded that the mean value on all of the orbits is also equal to φ 0 . However, in reality, on each orbit there might be two stars, one having φ = φ 0 + ∆ and the other one with φ 0 − ∆. Using only the map of the mean value one could not determine ∆. Having constrained the mean values of φ on the orbits, we can then project these quantities as a function of more accessible variables, such as spatial coordinates or velocities. In particular, in this contribution we are focusing on profiles as a function of position, and deprojected maps, but we could also plot φ as a function of velocity coordinates or actions.
Thus, our method, or any similar method using the same type of data, naturally produces a distribution function of φ at a given position which reflects the distribution of then mean φ on the orbits that pass through this position but it does not necessarily reflect the true distribution of φ there. To recover it, one should either add some additional, physically-motivated assumptions (like Zhu et al. 2020) or more data constraints.
Further natural question could be asked about possible degeneracies between stellar orbits. Our method is in fact based on the projected surface density of the orbits of each particle. If a projection of an orbit could be linearly decomposed into different orbits, then our technique would not be able to unambiguously assign stellar populations to it. In principle, one expects such an occurrence, since the orbital phase space can be labelled by three actions (see e.g. Binney & Tremaine 2008), while our data is inherently two-dimensional.
It is relatively easy to understand the issue for axisymmetric disc galaxies. In that case the orbits can be classified based on three conserved actions. The angular momentum L (equivalent to the azimuthal action) describes the size of the orbit. The radial action J r refers to the elongation of the orbit. The vertical action J z describes the vertical thickness of the orbit. First, let us consider only planar orbits (i.e. J z = 0). Then the projected density distribution of an orbit with a non-zero J r can be constructed as a sum of density distributions of circular orbits. On the other hand, the projected shapes of the orbits with a non-zero J z are clearly different, since they are oriented along the minor axis of the projection. We could concluded that the orbital projections are degenerate on the (L, J r ) plane. However, we are imposing physical limits on the values of the parameter φ i on a given orbit. Hence, we believe that the degeneracy is partially lifted since the model cannot freely assign the stellar population parameters.
The case of barred (i.e. non-axisymmetric) galaxies is more complicated. On one hand one still expects some degeneracy because of the relation between dimensions of the phase-space and the data. On the other hand, neither the classical planar orbital families x 1 -x 4 (Contopoulos & Papayannopoulos 1980), nor the vertically extended orbital families (Skokos et al. 2002), nor the orbits from the actual N-body simulations (e.g. Valluri et al. 2016;Gajda et al. 2016) exhibit obviously degenerate orbital projection. This issue certainly warrants further investigation.
Practical conclusions from our mock test, as well as from Zhu et al. (2020), are that there is plenty of information one can derive from spatially resolved mean maps of stellar population parameters. We showed that the radial and azimuthal profiles are very well constrained and one can retrieve also the vertical profile while applying some extra care.

M31 dynamical model and parameters used in stellar population modelling
As the basis for our stellar population fitting we employed the JR804 N-body model of the Andromeda Galaxy constructed by Blaña Díaz et al. (2018, see Tab. 1). It consists of 2 × 10 6 dark matter particles with Einasto density profile, 10 6 disc particles (which include the bar and its box/peanut bulge) and 10 6 classical bulge particles. Using M2M technique it was fitted to the moments of the velocity distribution derived by Opitsch et al. (2018) and to the 3.6 µm-band surface brightness maps from the Spitzer Space Telescope obtained by Barmby et al. (2006).
Blaña Díaz et al. (2017Díaz et al. ( , 2018 concluded that M31 is a barred galaxy. They found that the bar position angle (in the plane of its disc) with respect to the line of nodes is equal to 54.7 • . The bar rotates with a pattern speed of Ω p = 40 ± 5 km s −1 kpc −1 and has length of ≈ 4 kpc. The dark matter halo was found to follow an Einasto profile with a mass of 1.2 +0.2 −0.4 × 10 10 M ⊙ within 3.2 kpc and a stellar mass-to-light ratio in the 3.6 µm band of Υ 3.6µm = 0.72 ± 0.02 M ⊙ L −1 ⊙ , assumed to be a single constant. Blaña Díaz et al. (2018) in their paper provide a set of models that fit the data well and which could be thought of as "1σ models". Here we did not take those into account and thus our uncertainties does not include uncertainty of the dynamical model. Having a detailed dynamical model is an excellent basis for the modelling in this paper since the stellar orbits are thereby determined consistently with the photometry and kinematics, however some uncertainties remain. We assume that the population labels are mass-weighted and, hence, their best-fit distribution depend on the density distribution. Blaña Díaz et al. (2018) fitted the 3.6 µm IRAC image that traces the old giant stars (the bulk of the population). While the fitting was luminosity-weighted, the assumption of using a single constant mass-to-light ratio implies a trivial relation between particle masses and their 3.6 µm luminosity weights, which, in turn, implies that both mass-and light weighted averages give exactly the same result.
However, the measured kinematics are light-weighted in the V band where younger or more metal-poor stars could have a bigger impact in some regions (see e.g. Portaluri et al. 2017). Thus, while the surface mass distribution is well-constrained by the infrared photometry, the distribution of the particle orbits is biased towards the V band kinematics without considering the V band photometry. This could be a significant effect in the disk regions with relatively recent star formation. A future M2M model might be improved by fitting the V-band simultaneouly with the IRAC 3.6 µm to better model the kinematic structure in these regions.
An observant reader would notice also a small inconsistency in our approach. The Blaña Díaz et al. (2018) model assumes a constant M/L, however our approach posits a distribution of metallicities and α-enhancements. A reason not to vary M/L was our uncertainty whether the stellar population labels are precise enough to constrain the dynamics. Additionally, the projected variability of [Z/H] is of the order of 0.2 dex, which implies a difference in 3.6µm-band M/L of about 5% (Meidt et al. 2014) and about 10% in V-band (Vazdekis et al. 2010), both of which are of the same order as the ∼ 4% uncertainty of our model's M/L 3.6µm (including a systematic uncertainty from different choices of the dark matter profile). We conclude that the impact of a variable M/L is rather small and does not influence our main results. As a side note, we remark that to date none of the dynamical models of the stellar populations is proven to be fully internally self-consistent. Concerning the model of Poci et al. (2019), the remaining question is whether the mass distribution inferred from the orbital light-weights and their respective mass-to-light ratios reproduces the deprojected stellar mass distribution used to generate the orbits. Whereas  Zhu et al. 2020; "mass-to-clump ratio" in Portail et al. 2017a,b). It appears that making an actual fully self-consistent model is a worthwhile goal for a future research.
Following Blaña Díaz et al. (2018), we set 1 iteration to 1.18 × 10 −4 Gyr. Similarly, we fix the smoothing time scale to τ = 1.6 × 10 3 it, which corresponds to the orbital timescale at R = 5 kpc. We tested different values of ǫ that controls the strength of the force-of-change in (4). We found that in the case of [Z/H] the best results (i.e. the lowest χ 2 ) are obtained for log 10 ǫ = −1.9 and in the case of [α/Fe] for log 10 ǫ = −2.1. We smooth the observables initially for 3 × 10 3 it, we fit for 50 × 10 3 it, until χ 2 converges to a constant value and then we let the system relax for 7 × 10 3 it. As usual in the M2M modelling, after we finish fitting χ 2 increases in the relaxation phase, by about ∆χ 2 red ∼ 0.1, and stabilizes at the new and final value.

M31 metallicity
First we present our modelling of the [Z/H] distribution in Andromeda. In Figure 4 we show the data (top) and the fitted model (bottom). The white spaces correspond to the lines-of-sight not covered by the data. The (R x , R y ) reference frame (the same as in A&A proofs: manuscript no. unravelling_m31 Blaña Díaz et al. 2018) is rotated by 50 • clockwise with respect to the sky coordinates, so that the projected major axis of the M31 disc is almost aligned with the R x -axis. The kpc labels correspond to the sizes on the sky at the M31 distance. To convert R y into a distance in the plane of the disc, one should multiply it by (cos i) −1 ≈ 4.4. We also note that the data covers nearly the entire bar length. The possible values of [Z/H] are limited to the range -2.25 to 0.67 dex, as in Saglia et al. (2018). Overall, the fit seems to be very good (χ 2 red = 0.94), especially in the central part. However it is noticeably worse in the R y > 0 ′′ part, where we seem to underestimate [Z/H]. The data itself exhibit a strong asymmetry, with the top part of M31 having significantly higher metallicity than the bottom part (see top panel of Fig. 4). We hypothesise that this might be related to dust obscuration, since the R y > 0 ′′ part is closer to us. It is possible that the bulge stars are obscured by the star-forming disc, hence the ages are biased towards the younger ones, which in turn, through the well-known age-metallicity degeneracy, leads to an overestimated metallicity for the entire line-of-sight. On the other hand, we cannot exclude that the outer disk of M31 is highly metal-rich, despite the results of Gregersen et al. (2015). For a further discussion, see Saglia et al. (2018).
As described in Sec. 3, we have tried different values of G and N (see eq. 5). We found that the best (in terms of the final χ 2 ) parameters of the initial vertical profile are G = −0.14 dex/kpc and N = −0.01 dex. We stress again that these are not directly related to the final vertical profile.
In Figure 5 we present the M31 metallicity profiles, calculated from the fitted model. We plot [Z/H] as a function of cylindrical radius R (i.e. measured in the plane of the disc), height above the disc plane z and azimuthal angle ϕ. We took into account only the particles with |z| < 3 kpc. For the vertical profile we restricted R < 5 kpc and for the azimuthal profile we took into account only the particles with R < 5 kpc and |z| < 3 kpc. The uncertainty bands include both the uncertainty stemming from re-fitting the model to different realisations of the data within their errors and from using different vertical profiles within their respective 1σ uncertainties.
The first thing one notices in the radial profile is a spike in the central ∼ 200 pc, which is also discernible in the raw data of Saglia et al. (2018, Fig. 22) and was interpreted there as a sign of the classical bulge. Then, [Z/H] decreases further, up to ≈ 2 kpc, which later on we interpret in terms of the metallicity desert perpendicular to the bar. Further out, the azimuthally averaged profile starts to increase again, reaching a maximum around 4 kpc from the center. This coincides with the size of the bar in the Blaña Díaz et al. (2018) model. Still further out, in the disc outside the bar, the profile decreases once again, but does not exhibit a steady negative gradient, probably due to insufficient data or dust impact (Fig. 4). We decided to limit the radial range of the profile to 7 kpc, because of the asymmetry of the data (see Fig. 4), which starts to influence the results excessively beyond that point. Our profile can be compared to the radial metallicity profile obtained by the PHAT survey (Gregersen et al. 2015). Their radial profile has a metallicity maximum at R ≈ 4.5 kpc and further away [M/H] decreases linearly with radius. While we do not find an obvious linear decrease beyond the end of the bar, our mean values are consistent with the median values reported in their Fig. 9.
The azimuthal profile supports the Saglia et al. (2018) conclusion that [Z/H] is enhanced along the bar (which is aligned here with ϕ = 0, ±π). The Gregersen et al. (2015) [M/H] map also indicate that metallicity is enhanced along the bar, particularly close to its tip. The average vertical profile is con- sistent with linearity. Thus, we fitted a linear function using uncertainty-weighted least squares and obtained a vertical gradient of ∇ z [Z/H] = −0.132 ± 0.006 dex/kpc (see Fig. 5). Note that the uncertainty of the slope is likely underestimated, since it is driven by small uncertainty at |z| < 1 kpc. In Figure 6 we show deprojected maps of the mass-weighted mean metallicity. The minor axis of the disc is aligned with the z-axis, while the bar major axis is here aligned with the x-axis.
To create this maps we used only the particles within R < 7.2 kpc and |z| < 3 kpc. Then, they were time-smoothed using equation (2). To guide the eye, the mass density contours are overplotted in black.
In the central 3 kpc of the face-on view one can see an enhancement along the bar and depressions in the direction perpendicular to it. Further away one can notice maxima at both ends of the bar, around x ≈ 4−5 kpc. There is also a noticeable asymmetry between positive and negative y, which is a direct consequence of the top-bottom asymmetry of the data (Fig. 4).   We made an additional test to check whether this feature is dynamically stable, running the relaxation phase for 5 × 10 4 it, corresponding to ≈ 6 Gyr. We found that this feature is persistent, presumably because banana-shaped orbits around one of the Lagrange points L 4 /L 5 are assigned exceptionally high metallicity values. If the asymmetry is dust-related then one can view this feature as the model being too flexible. Our method treats the high [Z/H] at R y > 0 ′′ as physical, hence it assigns high Z to the model particles in the y > 0 region. Since the region at large, positive y may or may not be influenced by the dust, we believe that the bottom part (y < 0) is more trustworthy. It exhibits a noticeable ring of enhanced metallicity, beyond which [Z/H] decreases outwards, compatible with Gregersen et al. (2015). The edge-on view clearly exhibits an X-shape, caused by high-metallicity particles on vertically extended orbits in the boxy/peanut bulge. We note that the b/p bulge in the dynamical model extends to |x| ≈ 2.5 kpc. The X-shape persist further outwards, possibly driven by the metallicity increase in the plane up at x ≈ 4 kpc. It is possible that the metal-rich disk of M31 was thickened due to the recent merger (Bhattacharya et al. 2019). In the end-on view, the most prominent feature is the asymmetry with respect to the y = 0 kpc line, which we discussed above.

Enhancement of [α/Fe] in M31
We now focus on the [α/Fe] enhancement, which provides a signature of star-formation timescales (e.g. Tinsley 1979;Matteucci & Greggio 1986;Ferreras & Silk 2002;Thomas et al. 2005), i.e. larger α-enhancement signals quicker star formation. In Figure 7 we plot, as before, both the Saglia et al. (2018) data and our fitted model, which has an overall χ 2 red = 0.91. The possible values of the particle [α/Fe] are limited to the range from −0.3 to 0.5 dex, as in Saglia et al. (2018). As before, we tried a range of initial vertical profiles (see eq. 5). In this case the initial conditions with the lowest final χ 2 are given by G = 0.17 dex/kpc and N = 0.47 dex. At first glance we can see that the center of M31 is relatively more α-enriched and the shape of this feature is elongated along the projected minor axis of the galaxy. It is so because the lower-α bulge material is elongated along R x (see Saglia et al. 2018). Furthermore, [α/Fe] seems to decrease along the projected major axis for |R x | > 2 kpc, but is quite flat in the R y direction.
The [α/Fe] profiles, as a function of the cylindrical coordinates R, z and ϕ, calculated from our fitted N-body model are shown in Figure 8. We took into account only particles within |z| < 3 kpc. For the vertical profile we restricted R < 5 kpc, while for the azimuthal profile we considered only the region R < 5 kpc, |z| < 3 kpc.
In the inner 1.5 kpc the radial profile is rather flat, then it decreases, attaining a minimum at R ≈ 4 kpc. Further away, the profile gently bends upwards, while being consistent with a constant value. The drop in the innermost ∼ 300 pc can be attributed to presence of young stellar population . It is worth noticing that the data are consistent with the M31's disc having super-solar α-enhancement, and, hence, short star formation timescale (Thomas et al. 2005). Conversely to the metallicity, in the azimuthal profile α seems to be higher in the direction perpendicular to the bar.
We made also a linear fit to the vertical [α/Fe] profile and obtained a value of ∇ z [α/Fe] = (−0.007 ± 0.003) dex/kpc. Judging from the error band, [α/Fe] appears to be constant in the direc-A&A proofs: manuscript no. unravelling_m31 tion perpendicular to the disc, at least on average in the inner 5 kpc.
In Figure 9 we show deprojected maps of the mass-weighted mean α-enhancement. Similarly to Figure 6, the bar major axis is aligned with the x-axis and the disc rotation axis coincides with the z-axis. As before, the maps were time-smoothed following eq. (2) In the face-on view, the orientation of the [α/Fe] enrichment is clearly perpendicular to the bar direction and coincides with the metallicity deserts discussed in the previous subsection. Curiously, the face-on map exhibits a lower-α ring at R ∼ 3 −4 kpc, coincident with the metallicity enhanced ring visible in the bottom part of the left panel of Fig. 6. In the central parts of the edge-on and the end-on view the enhancement is high, with small-scale patterns which may or may not be real. In the outer part of the disc, the α-enhancement is smaller closer to the disc plane and grows with height. The tendencies in the central parts are most probably related to the classical bulge and the regions where it dominates over the disc population. The trend in the outer parts can be interpreted as an indication of the presence of a thick disc. The positive vertical gradient of [α/Fe] outside of the centermost part of the bulge is also reported for the Milky Way (e.g. Ness & Freeman 2016). Saglia et al. (2018) Several features of the stellar population distribution in M31 were clear from inspection of the maps obtained by Saglia et al. (2018), such as the [Z/H] enhancement along the bar and the [α/Fe] peak in the central 1 kpc. In addition, Saglia et al. (2018) constructed a simple decomposition of their stellar population data into various components (classical bulge, b/p bulge, bar and disk), which was based on using different regions of the sky to constrain the classical bulge and bar components in the mass model of Blaña Díaz et al. (2018).

Comparison to
It is instructive to compare our detailed results to the mean profiles of their model (blue and black lines in Figs. 21-24 of Saglia et al. 2018). We found a steep increase of metallicity and steep decrease of [α/Fe] in the inner 200 pc (50 ′′ ), which was directly mandated by the data we used. Saglia et al. (2018) ascribed the [Z/H] peak to the classical bulge, while the α decrease was attributed to the recent star formation. We agree with the Saglia et al. (2018) model that along the bar metallicity is constant, while [α/Fe] decreases. They found almost flat metallicity profiles for the b/p bulge and the disc . Our model showed that the b/p bulge has an X-shape in [Z/H], as one can find in Fig. 6. Thus, the metallicity profiles in the b/p bulge are growing with radius, up to its end at ≈ 2.5 kpc, at a range of heights z ∼ 0.3 − 1 kpc. The [Z/H] morphology of the disc shows diverse features: metallicity deserts perpendicular to the bar at |y| ∼ 2 kpc, then further away an enriched elongated ring at R ∼ 3 − 5 kpc and a possible further decrease beyond 5 kpc. The vertical distribution of [Z/H] beyond the b/p bulge seems to be flaring. Our [α/Fe] distribution in the central part of M31 show enhancement perpendicular to the bar and up to a significant height. Our model show a steeper α-overabundace decrease up to R ∼ 4 kpc than in Fig. 24 of Saglia et al. (2018).
It would be tempting to perform here a decomposition into classical bulge, b/p bulge, bar and disk components, using the model by Blaña Díaz et al. (2018) that explicitly consists of two types of particles (belonging to the classical bulge and disc). However, because of the spin-up through angular momentum transfer by the bar (Saha et al. 2012(Saha et al. , 2016, some of the classical bulge orbits overlap with those of the bar. Our current method cannot differentiate between the mean metallicities of the two components on the same orbit. From our final model, we can only tentatively say that the bulge component has a very steep gradient, while the profile of the bar/disc component is flat or even decreasing towards the centre, however we are not able to give uncertainties. Improving this aspect of the model seems possible but is beyond the scope of this paper.
Compared to the earlier modelling, several new qualitative results have emerged from our full particle modelling of the stellar population data: the discovery of a metallicity X-shape of the b/p bulge and of an [Z/H]-enhanced, low-α ring, as well as metallicity maxima at the tips of the bar. Furthermore, we found indications of an α-enhanced thick disc in M31.
Article number, page 10 of 14 Gajda et al.: Unravelling stellar populations in the Andromeda Galaxy

Comparison with other galaxies
The impact of the bar on the stellar population gradients has been an active topic for some time now. It has been found that the profiles along the bars are flatter than those perpendicular to bars Fraser-McKelvie et al. 2019;Neumann et al. 2020). Interestingly, the results of Fraser-McKelvie et al. (2019) indicate that gradients may be either positive or negative (see also Pérez et al. 2009), but the signs of the gradients in both directions remain correlated. However, judging from our M31 results, it is important to go beyond simple linear fitting, since the profiles in the bar regions can be more complicated. Looking at Fig. 5, one would infer a negative radial gradient in the central ∼ 2 kpc and a positive one further out, up to the bar end. It appears that the metallicity enhancement close to the bar ends was not yet widely appreciated. However, a closer inspection of available spatially resolved data sometimes reveals such trends. Many of the BaLROG bars show this feature (Seidel et al. 2016), as well as some of the galaxies analysed in the TIMER project (Gadotti et al. 2019;Neumann et al. 2020). Moreover, the Bovy et al. (2019) maps of the metallicity in the Milky Way also can be interpreted as if [Fe/H] is either enhanced at the bar end or in a nearby ring. On the other hand, Wegg et al. (2019) found more enhanced metallicities in their bar fields at 4 kpc central distance.
It is also interesting to compare our side-on view to other published examples. The edge-on galaxy NGC 1381 of Pinna et al. (2019b) shows clear signs of boxy/peanut shape both in metallicity and [Mg/Fe]. Williams et al. (2011) suggested that this galaxy may have a small classical bulge, which would make it different from M31 with its sizeable classical bulge. Pinna et al. (2019a) studied another galaxy with a b/p bulge, FCC 177. However, this one does not have obvious signs of the b/p bulge in metallicity and certainly not in [Mg/Fe]. However, curiously, [Fe/H] seems to have maxima in the disc plane at ≈ 3 kpc from the galaxy center, beyond the b/p bulge (of size ∼ 1.5 − 2 kpc, judging from the isophotes). In case of the Milky Way the model of Portail et al. (2017b, Fig. 10) suggests that [Fe/H] distribution in the centre has a peanut-like shape. The APOGEE data indicate presence of a metallicity maximum close to the bar end at l ≈ 30 • (e.g. Ness & Freeman 2016).
In the Milky Way α-enhancement is anti-correlated with metallicity, broadly speaking (e.g. Hayden et al. 2015, and many others). In our model this is also true to a certain degree. In the face-on view α-rich regions coincide with the metallicity deserts perpendicular to the bar. The high-Z ring is relatively less α-enhanced. In the edge-on views, at large distances from the center, a negative metallicity gradient correspond to a positive α gradient. Apparently, the center of M31 does not follow this trend, presumably due to the metal-rich, high-α classical bulge.

Origins of the spatial trends
Many of the models of the stellar populations origin in disk galaxies are focussed on the Milky Way, because it is the galaxy for which we have the most detailed data. However, these models are often useful for M31 as well, since it is a disk galaxy of similar mass and also includes a bar and a b/p bulge component. Often, the authors assume a certain initial metallicity distribution and test the impact of the further secular evolution (e.g. Martel et al. 2013;Martinez-Valpuesta & Gerhard 2013;Di Matteo et al. 2013). The other approach is tracing stellar populations in hydrodynamical simulations with star formation and feedback (e.g. Grand et al. 2016;Debattista et al. 2017Debattista et al. , 2019aTissera et al. 2016Tissera et al. , 2017Tissera et al. , 2019Taylor & Kobayashi 2017).
Di Matteo et al. (2013) analysed a case where a disc galaxy has initially an axisymmetric distribution of mass and metallicity. Initially, they assigned [Fe/H] such that it decreased linearly with the radius. During the further evolution a bar and spiral arms formed. The creation of the bar induced outwards radial migration and effectively mixed the stellar metallicities from the galaxy center up to the bar length, creating an enhancement along the bar and metallicity-deserts perpendicular to it (see also Debattista et al. 2019b, Fig. 15). This resulted in an azimuthal variation of metallicity. They found that the ratio of the azimuthal variation δ [Fe/H] (measured e.g. in the bar region, see Di Matteo et al. 2013, Fig. 8 and Sec. 3.3) and the initial metallicity gradient ∆ [Fe/H] can be approximated by δ [Fe/H] /∆ [Fe/H] ∼ 1 ± 0.5, independently of the initial gradient. In our case δ [Z/H] ≈ 0.02 dex (Fig. 5), which would correspond to an initial radial gradient in Andromeda of the order of ∇ R [Z/H] ∼ 0.02 ± 0.01 dex/kpc, compatible with the metallicity gradient in the outer disc of M31 derived from the PHAT survey (Gregersen et al. 2015), as well as typical radial gradients in observed galaxies (Sánchez-Blázquez et al. 2014) and in hydrodynamical simulations for redshift z < 1 (Tissera et al. 2017). Note that the exact shape of the bar-related enhancement depends on the initial mix of the stellar populations (Khoperskov et al. 2018). Further work is needed to understand the influence of the purported recent merger on M31 (Hammer et al. 2018;Bhattacharya et al. 2019) and whether it deposited enough material in the bar region to significantly change the mean stellar population parameters there.
Alternatively, the enhancement along the bar may be caused by the so-called kinematic fractionation effect Fragkoudi et al. 2017), by which the population that is colder before the bar formation would become more strongly aligned with the bar afterwards. If the colder population was additionally more metal-rich, this would result in a metallicity enhanced bar. In the Milky Way the more metal-rich populations exhibit distributions that are more strongly barred (Portail et al. 2017b). Moreover, the bar metallicity can be further enhanced during its growth, if the bar can capture new stars that are more metal rich (e.g. Aumer & Schönrich 2015).
In the case of the α-enhancement in M31 we did not find any relation with the bar. The high value implies that the stars in the entire bar region must have formed rapidly at early times (Pipino et al. 2006(Pipino et al. , 2008. The stellar population that formed the bar also appears to not have had an initial [α/Fe] radial gradient, i.e. the entire bar and the surrounding disk stars must have formed rapidly at early times. The original galaxy also should not have had multiple discs of different α, such as a distinct αrich thick disc, since this would probably also result in a barrelated non-axisymmetry (Fragkoudi et al. 2018).
The X-shaped side-on view of the metallicity is an expected outcome of the hydrodynamical simulations, both in the isolation ) and in the cosmological context (Debattista et al. 2019a;Fragkoudi et al. 2020). The X-shape in the region of our model's b/p bulge (|x| 2.5 kpc) is probably somewhat weaker than in the aforementioned works, since the dynamical b/p bulge is weaker due to the presence of the classical bulge. Differently than in the aforementioned models, the metallicity distribution in M31 keeps flaring with growing distance from the center. We are not entirely sure what is its origin of this effect and how significant it is. Certainly, more spatially extended data beyond the bar end would help to constrain the model more strongly. The side on-view of the [α/Fe] in M31 looks different from what was obtained in Auriga simulations . Note that these Auriga galaxies do not have classical bulges (as inferred from small Sersic indices), while M31 clearly has one. Note also that their average [α/Fe] are much closer to the solar value than in Andromeda.
Metallicity enhancement along bars and maxima at their ends are also found in hydrodynamical simulations (Debattista et al. 2019a;Fragkoudi et al. 2020). In particular, two of the runs (Au18 and Au23) analysed by Fragkoudi et al. (2020) exhibit both [Fe/H] maxima at the bar ends and rings of enhanced metallicity and relatively lower-α. The authors of that study relate the rings to on-going star formation in the region between the bar end and the corotation radius. The runs Au18 and Au23 are different with respect to the others analysed by Fragkoudi et al. (2020) also in terms of their rotation curves, which are not centrally peaked, similarly to our model of M31 (Blaña Díaz et al. 2018) and the Milky Way (e.g. Portail et al. 2017a). Those bars are also significantly smaller than their corotation radii, i.e., they are relatively slower. The enhancement close to the tip of the bar could be related to the star formation activity observed there in some of the galaxies (e.g. Phillips 1996). Debattista et al. (2019b, Fig. 15) found that if one tags the particles with metallicity corresponding to an axisymmetric galaxy and let the galaxy form a bar, then the metallicity becomes enhanced along it, but the maxima at the bar ends do not form.
High metallicity rings and arcs beyond the metallicity desert have also been discerned in Khoperskov et al. (2018), where they are related to spiral arms (see also Debattista et al. 2019a, Fig. 8). They arise from stellar populations with different kinematics and different metallicity patterns, which respond differently to the spiral arms. The ring we find in our model in principle could have a similar origin.

Conclusions
We devised a new chemodynamical technique based on the made-to-measure (M2M) modelling framework, enabling us to constrain the distribution of the galactic stellar populations. As an input, the method uses mass-weighted maps of a mean stellar population parameters, e.g. metallicity, age or α-enhancement. We start with a dynamical N-body model, constrained by surface density and kinematics. Then, we tag the particles with a single value of e.g. metallicity, corresponding to mean metallicity along the orbit, and then we adjust those values to match the observational constraints. We tested our method on mock data and found that the profiles in the plane of the galaxy (i.e. radial and azimuthal) are well reproduced, while in the case of the vertical profile we are able only to obtain an average gradient, without finer details.
We applied our technique to the Andromeda Galaxy. We used [Z/H] and [α/Fe] maps derived by Saglia et al. (2018) to constrain the three-dimensional distribution of metallicity and α-enrichment in the context of the dynamical model of Blaña Díaz et al. (2018) for M31.
Our main results can be summarized as follows: 1. We find that the metallicity is enhanced along the bar, while perpendicular to it we see a clear depression, akin to the so-called star formation deserts (James et al. 2009;Donohoe-Keyes et al. 2019). We find also that the metallicity distribution exhibits a peak at the radius corresponding to the bar size. In the face-on view we can discern an elongated, [Z/H] enhanced ring. 2. Closer inspection of the side-one view reveals that the [Z/H] enhancement has an X-shape caused by the boxy/peanut bulge. We estimate the average vertical gradient of [Z/H] to be −0.132 ± 0.006 dex/kpc. While typical values for the vertical gradient are about twice larger than what we find for M31 (Molaeinezhad et al. 2017), for example in NGC 1381 (Williams et al. 2011;Pinna et al. 2019b) and in the MW (Ness et al. 2013;Rojas-Arriagada et al. 2014), there are some galaxies with similar vertical metallicity gradients in the sample of Molaeinezhad et al. (2017). 3. On average, the [α/Fe] enhancement in Andromeda is high compared to the solar value, indicating short formation timescales. The centre of M31 is clearly more α-enhanced than the surrounding disc, probably due to the presence of the classical bulge. We also find a relatively lower-α ring, corresponding to the metallicity ring. The average vertical gradient is flat, however a closer look reveals more structure.
In the centre α-enhancement decreases slightly, while in the outer parts it increases with height.
The central behaviour is an imprint of the classical bulge, while the outer distribution may signify a presence of a more strongly α-enhanced thick disc. From the data of Saglia et al. (2018) and our model the following formation pathway for the Andromeda Galaxy can be deduced. Given the overall high levels of α-enhancement, both the bulge region and the inner disc must have formed relatively quickly, with faster star-formation in the bulge region. Since the current structure of the metallicity and α distribution are qualitatively different, they must have been different initially. The distribution of [α/Fe] is consistent with no radial gradient or any other structure in the original disc. The galaxy assembly resulted in a sharp peak of metallicity in the central few hundred parsecs and a more gentle negative gradient in the remaining disc. The formation of the bar lead to a re-arrangement of the [Z/H] distribution, causing the formation of metallicity-deserts and a flat gradient along the bar. Afterwards, the star formation continued close to the bar ends in the leading edges of the bar, producing metallicity enhancements in the ansae and the [Z/H] enhanced, lower-α ring. The star formation in the very center at ∼ 200 pc also continued, however it was rather mild, leading to only a small decrease of [α/Fe] there. Andromeda probably experienced recently a fairly massive minor merger (Hammer et al. 2018;Bhattacharya et al. 2019), which most probably lead to flattening of the gradients (e.g. Taylor & Kobayashi 2017). It would be interesting to further investigate the impact of a merger on the distribution of the stellar populations in a barred galaxy.
In this work we constrained our models with spatially resolved maps of line-of-sight averages of a given stellar population parameter. Thus we could only derive the mean metallicity and [α/Fe] in the corresponding part of the phase-space (deprojected space). The obvious way forward is to use metallicitydistribution-resolved maps, which would enable us to study various properties of the galaxies as function of metallicity.
Our newly developed technique proved to be a valuable tool for studying the distribution of the stellar populations in galaxies. To date, analysis was always confined to looking at galaxies either face-on or edge-on. Here we are able to model the threedimensional distribution and thereby connect the results from both perspectives. In the future, we plan to apply our technique to other galaxies with spatially-resolved maps of metallicity, αenhancement and age, thereby obtaining a more complete picture of the chemodynamical structures in disk galaxies.