THE LUMINOSITY PROFILE AND STRUCTURAL PARAMETERS OF THE ANDROMEDA GALAXY

The library of SB profiles for M31 is vast. Table 1 presents a list of the major photographic and digital databases for luminosity profiles and star counts of M31. The SB profiles extracted along the major axis from the galaxy images listed in Table 1 are shown in Figure 1 (star counts are presented in Figure 5).

Zoom In Zoom Out Reset image size

Table 1. Sources of Surface Brightness Profiles for M31

Reference	Bands	Pixel Scale	Minimum Experiment Time
Walterbos & Kennicutt (1987)	UBVR	Photographic plate	300–7200 s
Pritchet & van den Bergh (1994)	V	Star counts	2700 s
Choi et al. (2002)	I	203 pixel−1	2400 s
Irwin et al. (2005)	Gunn i	Star counts	800–1000 s
Beaton et al. (2007) 2MASS 6X	JHK	100 pixel−1	46.8 s
Barmby et al. (2006) IRAC	3.6 μm, 4.5 μm	0863 pixel−1	62–107 s
Gilbert et al. (2009)	V	Spectroscopic star counts	3600 s

Download table as:  ASCIITypeset image

The much-cited library of UBVR luminosity profiles of M31 from photographic plates by Walterbos & Kennicutt (1987) already highlights the so-called "10 kpc arm," a noticeable bump above a rather exponential disk profile. However, these optical light profiles do not sample the galaxy center well; the starlight in the UBV filters is also sensitive to dust extinction and mostly dominated by bright massive stars, which contribute only a small fraction of the total galaxy stellar mass. KB99 combined and analyzed various V-band light profiles of M31. However, we focus here on redder, less extinction-prone and more spatially extended luminosity profiles.

Choi et al. (2002, hereafter Choi02) derived 17 × 5° B- and I-band CCD mosaics of M31, including M32 and NGC 205, with the Kitt Peak National Observatory Burrell Schmidt telescope. Each mosaic included 35 fields, with a typical exposure time of 40 minutes per field in I. An azimuthally averaged SB profile of this I-band mosaic was extracted via isophotal fitting by Worthey et al. (2005); see the red line in Figure 1. For the purpose of the current analysis, we used the same I-band mosaic from Choi02 to extract a revised I-band azimuthally averaged SB profile as well as I-band major- and minor-axis wedge cuts (see Section 2.2). We reassessed the absolute calibration of the Choi02 data by performing aperture photometry of field stars in that mosaic and adjusting the stellar magnitudes to those computed by the Sloan Digital Sky Survey (SDSS). We used the same point-spread function fitting routines for the Choi02 and SDSS star fields. This procedure yielded a photometric accuracy of ∼0.1 mag (McDonald et al. 2009, hereafter McD09). The calibration from SDSS stars led to a readjustment of 0.4 mag from the photometry reported in Worthey et al. (2005).

Figure 1 shows that the Choi02 data extend as far as the Walterbos & Kennicutt (1987) profiles and better sample the central regions. For this reason, and because the rest-frame I band is a better tracer of the stellar mass than UBV, we use the Choi02 data as our principal anchor for the derivation of a composite light profile for M31. While it would be desirable to combine the data presented in Table 1 into one global composite profile, color gradients thwart that project. Different bands yield different structural parameters, owing to differently distributed stellar populations and dust, selective dust extinction, etc.

High-resolution, wide-field infrared mosaics of M31 remain scant due to the smaller size and lower efficiency of IR arrays, as compared to optical CCDs, and the significant complications due to the rapidly varying infrared sky background. Existing infrared maps come from the Two Micron All Sky Survey (2MASS) 6X survey (Beaton et al. 2007) and the Spitzer Infrared Array Camera (IRAC) survey (Barmby et al. 2006). Beaton et al. (2007) produced a 2.8 deg² JHK map of M31 as part of the 2MASS 6X program.⁸ The exposure time per source is six times that of the nominal 2MASS integration, or six times 1.3 s per frame times six frames per source, that is, 46.8 s per source. We have extracted near-IR SB profiles from the 2MASS 6X JHK mosaics according to standard surface photometry techniques (e.g., Courteau 1996; McD09). Counts were calibrated to magnitudes using 2MASS background stars. SB levels were inferred from the 2MASS 6X images using the following transformation: SB[mag arcsec⁻²] = −2.5log (I[counts]) + z.p., where the zero points are J:22.81, H:22.35, and K:21.85 mag arcsec⁻². The H- and K'-band 2MASS 6X azimuthally averaged SB profiles are reported in Figure 1 as brown lines. The 2MASS data (both LGA and 6X) remain ill-constrained beyond the 10 kpc arm due to uncertain sky subtraction (Barmby et al. 2006). Indeed, while the proximity of M31 enables a sharper view of its stellar content, its large apparent size is detrimental for the reliable measurement of the sky level far from the galaxy on timescales shorter than those of the natural fluctuations of the IR sky. Producing an accurately sky-subtracted image for M31 is an extremely challenging task that we address in a forthcoming publication (J. Sick et al. 2011, in preparation). We therefore leave out the 2MASS data from our analysis.

The Spitzer mosaic of Barmby et al. (2006) covers 37 × 16 in the four IRAC bands at 3.6, 4.5, 5.8, and 8 μm. These mosaics sample the galaxy light with 0861 pixels from the center out to the "10 kpc" arm. The integration times per pixel ranged from 62 s in the disk to 107 s in the outskirts. Our extracted IRAC azimuthally averaged SB profiles are shown as gray-to-black lines in Figure 1. For our study, we discarded the 5.8 and 8 μm profiles which are spoiled by polycyclic aromatic hydrocarbon (PAH) emission, especially near the 10 kpc arm. The stellar light at 3.6 and 4.5 μm bands is significantly less contaminated by hot dust features. For simplicity, we only used the 3.6 μm band for our analysis; the 3.6 and 4.5 μm data yield similar profile decompositions.

Although the Spitzer/IRAC images do not have the superb resolution of HST, the IRAC luminosity profiles in Figure 1 already show the suggestive signature at R ≲ 15 pc of a nuclear component. We shall address the modeling of the nucleus in Section 3.1.

It is interesting to note that the signature of the 10 kpc arm shifts inward as a function of wavelength. This is expected if the blue stars lie at the front of the spiral density wave where they were formed; the dust (e.g., PAHs) that results from the evolution of those massive stars will naturally drift behind the stellar front, as observed clearly in the brightness enhancement between 5 and 9 kpc of the IRAC 5.8/8.0 μm light profiles in Figure 1.

The luminosity profiles shown in Figure 1 are all extracted from azimuthally averaged isophotal fits. However, as with most other spiral galaxies, the light of M31 is a superposition of multiple components with different ellipticities. Isophotal fitting software packages, originally designed for single-component elliptical systems,⁹ average out the ellipticity information of a galaxy at each pixel such that the fitted ellipses no longer represent one independent structure. This is perhaps the major drawback of profile decompositions based on azimuthally averaged luminosity profiles.

Fortunately, for the case of M31, the position angles of the bulge and disk are within 20° of each other, which lessens the mismatch of superimposed isophotes. However, the ellipticities of the two components are still quite different (see Section 7). Consequently, we shall also explore 1D profile decompositions based on the simultaneous modeling of minor- and major-axis wedge cuts. Such cuts then allow for the modeling of independent galaxian structures. We construct wedge cuts at different azimuthal angles as follows. In order to keep a high signal-to-noise ratio along the cut, we adopted square bins with variable width, W, given by

$$\begin{equation} W= N p (e^{n/N} -1), \end{equation} \tag{ 1 }$$

where n is the bin number, p is the minimum number of pixels at small radius, and N is a constant. The opening angle of the wedge is thus a nonlinear function of galactocentric radius. We have used p = 3 and N = 25.

Figure 2 shows the position of the major-axis wedge onto the IRAC 3.6 μm image. A median SB is calculated in each variable-width bin. To remove bright stars, all pixel fluxes more than 3σ above this value are excised and the SB profile is recomputed.

Zoom In Zoom Out Reset image size

The minor- and major-axis wedge cuts extracted from the Choi02 and IRAC 3.6 μm images are shown in Figure 3. The major-axis cuts are also compared with the azimuthally averaged SB profiles at the I band and at 3.6 μm. Despite our concerns that azimuthal averaging yields, in principle, an unrepresentative blend of distinct luminous components each with different ellipticities, the azimuthally averaged profiles for M31 differ only slightly from the precise cuts. The cuts and averaged profiles are thus interchangeable but the latter have higher signal-to-noise levels.

Zoom In Zoom Out Reset image size

Subtraction of the I-band and 3.6 μ cuts in Figure 3 reveals color gradients in the inner parts (R_maj < 20 kpc) of M31. These gradients, shown in Figure 4, are consistent with the bulge being generally populated by more metal-rich stars than the disk, as is seen in other Sb–Sc-like galaxies (MacArthur et al. 2004, 2009; Moorthy & Holtzman 2006; Roediger et al. 2011). The bluing of the disk, especially beyond the 10 kpc ring, confirms that the I-band and 3.6 μm luminosity profiles cannot simply be merged (without corrections) to build an extended composite light profile for M31.

Zoom In Zoom Out Reset image size

Rigorous error bars are crucial to the proper model decomposition of a luminosity profile. For azimuthally averaged SB profiles, the quoted intensity per radial increment (here, per pixel) is the median of all the pixel intensities along a given elliptical isophote (light contour) and the brightness error per radial bin is the standard deviation of all those pixel intensities along the isophote with respect to the median intensity value (Courteau 1996). For a wedge cut along any radius, the median intensity per bin is computed as in Section 2.2; the quoted error is the rms deviation about the median value of the pixels in each bin. The latter error estimates can never be an exact representation of the true galaxian (sky-subtracted) photon noise at each bin since each pixel element is smaller than the seeing disk, and all adjacent pixels are thus correlated. The size of the error bars per bin is, however, comparable to the amplitude of the point-by-point fluctuations (see Figure 3), suggesting that our error estimates are reasonable.

The error bars as quoted by the authors for the Irwin05 profiles are significantly larger, at large radii, than the point-by-point fluctuations. For our modeling of the Irwin05 data, we have recomputed those errors as the standard deviation about linear fits through neighboring data points. These errors closely reflect the point-by-point fluctuations of the data.

In addition to the galaxy images and surface brightness profiles, information about the shape of the galaxy components can be gleaned from star counts measured in fields around the disk of the galaxy where crowding effects are lessened. We have combined the extended star counts of the M31 stellar halo largely along the minor axis by PvdB94, Irwin05, and Gilbert et al. (2009a, hereafter G09). These star counts cover the range 20 kpc ⩽R_min ⩽ 150 kpc (along the minor axis; see also Ibata et al. 2007; McConnachie et al. 2009).

PvdB94 obtained digital star counts with the Canada–France–Hawaii Telescope in the V band for eight fields along M31's minor axis. Each field was observed three times for a minimum of 900 s per integration. These counts were then combined with photographic luminosity profiles by de Vaucouleurs (1958) and Walterbos & Kennicutt (1987) and digital imaging by Kent (1983). The PvdB94 minor-axis brightness profile is shown with red circles in Figure 5. Based on these data, PvdB94 postulated that the halo of M31 along the minor axis was well represented by a single de Vaucouleurs law over the range 0.2 kpc ⩽R_min ⩽ 20 kpc.

Zoom In Zoom Out Reset image size

Irwin05 used stellar counts from the Isaac Newton Telescope to determine M31's luminosity profile along its southeast minor axis. They combined PvdB94's data with faint red giant branch star counts to trace the minor-axis stellar distribution out to a projected radius of ∼55 kpc, where μ_V ∼ 32 mag arcsec⁻². They exposed typically 800–1000 s per field per bandpass in Johnson V and Gunn i. We show their Gunn-i minor-axis luminosity profile as green dots in Figure 5. Irwin05's analysis corroborated PvdB94 with a de Vaucouleurs profile out to a projected minor-axis radius of ∼20 kpc. Beyond that radius, Irwin05 surmised, the light profile would assume a more exponential shape with a scale length of 14 kpc. We will see in Section 6 that the M31 minor-axis cut is indeed a poor tracer of disk light. Indeed, our decompositions argue against a de Vaucouleurs profile (Sérsic n = 4) for the M31 bulge.

Guha05 presented a V-band SB profile along M31's southeastern minor axis based on counts of resolved member red giant stars in a Keck/DEIMOS spectroscopic survey. The spectroscopic exposure time with Keck/DEIMOS was 1 hr per field. Red giants in M31 were identified (and distinguished from foreground Milky Way dwarf star contaminants) using a combination of photometric and spectroscopic diagnostics (Guhathakurta et al. 2006; Gilbert et al. 2006). To estimate M31's stellar surface density as a function of radius, the observed ratio of M31 red giants to Milky Way dwarf stars was multiplied by the surface density of Milky Way dwarf stars predicted by the "Besançon" Galactic star-count model (e.g., Robin et al. 2003, 2004). This estimate was then calibrated to the V-band intensity by comparing it with the PvdB94 data in overlapping regions (also based on star counts). The minor-axis SB profile from V-band counts, presented here as black dots in Figure 5, includes two improvements relative to the Guha05 analysis—as reported in G09: (1) a larger number and wider spatial distribution of Keck/DEIMOS spectroscopic fields providing denser sampling over the radial range 10 < R < 165 kpc (Kalirai et al. 2006; Gilbert et al. 2007), and (2) correction for field-to-field variations in the spectroscopic sample selection function.

The extended minor-axis composite profile of M31 shown in Figure 5 is a combination of the minor-axis wedge cut from the image of Choi02 (Section 2.2) and the deep minor-axis star counts of Irwin05, PvdB94, and G09. The star counts were all normalized to our recalibrated I-band profile for Choi02. Altogether, our composite profile yields a detailed picture of M31's global I-band luminosity distribution for μ_I < 30 mag arcsec⁻². Figure 5 may be compared, for R_minor < 80 kpc, with Figure 31 of Ibata et al. (2007), using V − I ∼ 2.

Our modeling of the composite profile can either treat all the data points independently ("Unbinned") or use a uniformly sampled profile that has been averaged over all overlapping data ("Binned"). Both techniques are often found in the literature and we wish to test them here. A weighted fit should in principle yield similar solutions for both the unbinned and binned profiles since the former will have more data points per bin (in regions of data overlap) while the latter will have, on average, smaller error bars per bin. Reality is not always so simple.

In order to create a single, composite, "binned" minor-axis profile from the four above-mentioned profiles, we adopted a logarithmic binning as in Equation (1) with p = 3 and N = 25. These parameters were chosen to match the inner regions of the I-band profile minor-axis profile. The surface brightness per bin is the error-weighted average of the surface brightnesses from all of the available data in that bin. The error per bin is the root-mean-square of all the data point errors in that bin. We will analyze both the unbinned and binned composite profiles in Section 7.

As noted above, Lauer et al. (1993) were first to show with their HST/WFPC1 images that the M31 nucleus was double-peaked. Despite the poorer resolution of our IRAC images, and our inability to resolve two distinct sources, the gentle central rise of the IRAC light profiles for R < 15 pc in Figure 1 is clearly suggestive of a nucleus component. KB99 used Lauer's HST/WFPC1-2 images of the M31 core and high spatial-resolution CFHT/SIS spectra to infer that the M31 nucleus is a separate component with an origin different from that of the bulge. The well-resolved HST V-band images analyzed by KB99 show a nearly exponential nucleus (see Equation (4) below) with a scale length R_n = 088 = 3.3 pc (J. Kormendy 2011, private communication). The nucleus contributes only a very small fraction (<0.05%) of M31's total luminosity budget. M31's nucleus is indeed sufficiently small that the model parameters for the bulge, disk, and halo are largely unaffected by it. Therefore, while we embrace the results by Lauer et al. (1993) and KB99, we now ignore the nucleus component in our luminosity profile decompositions.

We assume that the 2D luminosity distribution of the M31 bulge has a projected profile given by the Sérsic (1968) profile,

$$\begin{equation} I_b(R)=I_e\exp \left\lbrace -b_n\left[\left({R\over {R_e}}\right)^{1/n}-1\right]\right\rbrace, \end{equation} \tag{ 2 }$$

where R_e is the projected half-light radius, I_e is the intensity at R_e, the exponent n is the Sérsic shape parameter, and b_n = 1.9992n − 0.3271 (Capaccioli 1989). With n = 1 or n = 4, the Sérsic function reduces to the exponential or de Vaucouleurs functions, respectively.

The total extrapolated luminosity of a Sérsic bulge is given by

$$\begin{equation} L_b=2\pi (b/a) \int _{0}^{\infty }I_{b}(R) R\,dR=2\pi (b/a) I_{e}R_{e}^{2}{{e^{b_{n} \,}n\Gamma (2n)}\over {b_{n}^{2n}}}, \end{equation} \tag{ 3 }$$

where a and b are the semimajor and semiminor axes of the bulge. We will use Equation (3) to compute the bulge total light fraction in Section 7. See Graham & Driver (2005) for more details about the Sérsic function.

We assume that the projected light profile of the disk is described by a basic exponential function,

$$\begin{equation} I_d(R)=I_0\,\exp \lbrace - {R / R_d } \rbrace, \end{equation} \tag{ 4 }$$

where I₀ is the disk central intensity and R_d is the projected disk scale length. The total extrapolated luminosity of an exponential disk is given by

$$\begin{equation} L_d=2\pi (b/a) I_{0} {R_d}^{2}, \end{equation} \tag{ 5 }$$

where a and b are now the semimajor and semiminor axes of the disk.

Various authors (Kormendy 1977; Baggett et al. 1998; Puglielli et al. 2010) have considered an inner disk truncation as an alternative to modeling galaxy light profiles with a dip at the bulge–disk transition. Such a dip can be understood through various effects, such as dust extinction (MacArthur et al. 2003, hereafter Mac03) or the mixing action of a dynamical hot stellar bulge or bar (Puglielli et al. 2010). However, we ignore the possibility of a core in the inner disk of M31, as the current study based on imaging data alone offers no leverage for testing this hypothesis. We limit our modeling of the disk to a pure exponential function with a constant scale length at all radii.

The deep star counts from Irwin05 and G09, shown in Figure 5, are clearly indicative of an additional component at R_min ≳ 12 kpc, which appears to be independent of the galaxy disk (see also Ibata et al. 2007; Tanaka et al. 2010). As shown by Gilbert et al. (2007, G09), the kinematics of the stars in tidal streams out to 60 kpc show a large spread in velocities, consistent with expectations for a halo population. In Section 7.3, we also suggest that the M31 bulge and halo components are likely independent on account of their different spatial distributions.

We model the (1D) halo with either a Sérsic function (Equation (2)) or a power law:

$$\begin{equation} I_{h}(R) = I_* \left\lbrace \frac{1 + \left(R_*/a_h\right)^2}{1 + \left(R/a_h\right)^2}\right\rbrace ^\alpha, \end{equation} \tag{ 6 }$$

where R_* is a turnover radius and I(R_*) = I_*. The magnitude formulation is simply μ_h = −2.5log I_h. Based on visual examination of Figure 5, we adopt R_* = 30 kpc (measured along the minor axis) where the halo is dominant. The value of R_* is arbitrary and only affects I_*. The quantity a_h sets the amplitude of the profile in the inner parts. Note that at large radii, the Hernquist model (Hernquist 1990) and the Navarro–Frenk–White (NFW) function (Navarro et al. 1996), as explored for the M31 halo by, e.g., Ibata et al. (2007) and Tanaka et al. (2010), are special cases of our power-law model.

For the power-law halo, the total luminosity is

$$\begin{eqnarray} L_h&=& { 2\pi (b/a) I_{*} {R_*}^{2}\frac{1 + a_*^2}{2(\alpha -1)}}\nonumber\\ && \times\, { \left\lbrace \left(\frac{1 + a_*^2}{a_*^2}\right)^{\alpha -1} - \left(\frac{1 + a_*^2}{s_{\rm max}^2 + a_*^2}\right)^{\alpha -1} \right\rbrace }, \end{eqnarray} \tag{ 7 }$$

where a_* = a_h/R_* and s_max = R_max/R_*. In this work, we integrate out to a maximum radius R_max = 200 kpc along the minor axis.

We also assume a spherical halo (a = b) as the current star counts offer no leverage on the halo shape. If M31's halo is anything like the Milky Way's, our assumption of a spherical halo is well justified (Majewski et al. 2003; De Propris et al. 2010; however see Bell et al. 2008).

We can exploit different modeling techniques through a combination of 1D NLLS and MCMC analysis with "Binned" or "Unbinned" cuts. The "unbinned" method comes from our consideration of all the data points in our four individual cuts (Choi02, Irwin02, PvdB94, and G05) treated individually. The "binned" approach consists of averaging the four data sets into one composite profile; the overlapping data at any given radius are weight-averaged according to their respective errors. The differences in model parameters for the "Binned" versus "Unbinned" data, as seen in the Model pair S/T, can also be rather significant. Unless a precise argument can be made for data binning (e.g., improvements of the signal-to-noise ratio), the latter is best avoided.

Figure 8 shows the composite profile decompositions with Models U (power law; left) and W (Sérsic; right). Both models look statistically comparable given the similar residual distributions. Based on global χ-square statistics for NLLS Models R and V, which rely on minor-axis data only, and visual impression of the data-model residuals (figure not shown), the power-law fit seems to be a marginally better description of the outer profile of M31.

However, the MCMC "MinMaj" Models U and W, which must also account simultaneously for the major-axis I-band data of Choi02, are intrinsically more rigorous than Models R and V. The probability distributions for the parameters of Models U and W are also very similar. And while the fit residuals in Figure 8 are essentially equivalent, we see that the disk–halo transition for Model U occurs at R_min ∼ 8 kpc while for Model W the transition occurs at R_min > 10 kpc. This is made even clearer when we plot, in Figure 9, the relative light fraction of each component for Models U and W. As we will see in Section 7.4, kinematic constraints for M31 red giants (Gilbert et al. 2007; G09) suggest a sub-dominant disk for R_min > 9 kpc and thus Model W could be disfavored on those grounds alone. Furthermore, with its high-n index for a Sérsic halo, Model W calls for an excess of "halo" stars over disk stars in the inner parts of M31; the latter has yet to be detected (e.g., Saglia et al. 2010) though the scales involved are very small and both the measurements and their interpretation would be very challenging. For these reasons, we favor the power-law description (Model U) for the halo component.

The bulge of M31, long considered a prototypical de Vaucouleurs spheroid (e.g., Walterbos & Kennicutt 1988; PvdB94; Irwin05), is here confirmed to be less concentrated with a mean Sérsic shape index n ≃ 2.2 ± 0.3 for all the models and data reported in Tables 2–5.

Based on the IRAC luminosity profile at 3.6 μm, the mean bulge effective radius, R_e, projected along the major axis, is 1.0 ± 0.2 kpc, and the mean bulge effective SB, μ_e, is 16 ± 0.2 mag arcsec⁻². The surface brightnesses are, again, different for the I and IRAC bands.

Given the known covariances between the parameters R_e and μ_e for spheroid systems (Mac03), it is reassuring that the variations of three bulge parameters for the different NLLS and MCMC codes are so small. For small bulges, the Sérsic parameter n is not significantly covariant with other bulge parameters (Fisher & Drory 2010). Indeed, we find that n is independent of the decomposition method, data projection (minor versus major cuts), but it does depend slightly on bandpass (n being slightly larger at redder wavelength).

Our results are consistent with the earlier report by KB99 that n = 2.19 ± 0.10 and R_e = 1.0 ± 0.2 kpc for the M31 bulge (shown as the gray annulus in Figure 7). To our knowledge, KB99 is the first reported claim that the M31 bulge is less cuspy than a de Vaucouleurs profile (see their Section 2). Their V-band decomposition included a nucleus component but its impact on the bulge parameters, as we saw, is negligible.

Unaware of KB99, Worthey et al. (2005) decomposed the Choi02 I-band azimuthally averaged surface brightness profile using the Mac03 NLLS algorithm to find n = 2.0 ± 0.2 and R_e = 1.0 ± 0.3 kpc for the M31 bulge. This solution matches our Model N.

Seigar et al. (2008) used the Spitzer/IRAC 3.6 μm light profile and an NLLS model decomposition to find for the bulge n = 1.71 ± 0.11 and R_e = 1.93 ± 0.10 kpc. Their result can be compared with our Model F; see also the brown dot in Figure 7. They find an extended, brighter, dominant bulge with B/D = 57% ± 2%, whereas ours is sub-dominant with B/D ∼ 43%. We offer no explanation for their derivation of an unusually large bulge for M31.

As seen in Table 2, the bulge ellipticities determined by MCMC, with simultaneous major- and minor-axis analysis, and isophotal contours average _b = 0.25 ± 0.4 are clearly at odds with _b = 0.37 from GALFIT. This discrepancy is because, unlike our 1D profile analyses, GALFIT accounts for P.A. variations, and the inner isophotal contours are the sum of individual components (bulge and disk) with different ellipticities and P.A.s. The GALFIT P.A.'s for the bulge and disk differ by 20°, whereas our MCMC code assumes co-aligned B/D components. In Model Q, we see that forcing the MCMC ellipticities into GALFIT yields parameters that are closer to Model F. One can then only assume that a Markov Chain 2D decomposition (to be tested elsewhere) would recover results similar to Model E. Finally, a comparison of Tables 2 and 4 (Models L versus R and M versus U) shows that inclusion of the halo component does not affect the bulge light significantly (if at all, within the decomposition errors).

The spread of disk solutions from model-to-model variations might be expected to be larger than that of the bulge, in light of the inherent sensitivity from code-to-code to enhanced non-axisymmetric features. However, uncertainties in the disk measurements translate directly into bulge parameters. Typically, published disk scale lengths differ by ∼20% from author to author (Mac03). Figure 7, however, shows that our disk scale length solutions differ by less than that (at the 1σ level).

Our basic exponential disk model applied to the IRAC data has a scale length, R_d, of 5.3 ± 0.5 kpc. The rather small 10% uncertainty reflects the range of our model solutions from MCMC to 1D NLLS and 2D NLLS (GALFIT), i.e., from short to longer values. Proper accounting of model-data discrepancies, as with MCMC, yields shorter disk scale lengths. A realistic disk scale length for M31, measured in a band dominated by the older, low-mass disk stars, probably lies closer to 5 kpc.

For their B/D decomposition of V-band luminosity profiles, KB99 obtained R_d = 1698'' = 6.5 kpc. We know from Mac03 that V-band disk scale lengths are roughly 10% longer than I-band and 25% longer than H-band disk scale lengths. Thus, barring any additional color effect from H and 3.6 μm, the 6.5 kpc V-band scale length is equivalent to 5.2 kpc at 3.6 μm and 5.9 kpc at the I band in very good agreement with Models F and N. Worthey et al. (2005) measured an I-band value of R_d = 5.8  ±  0.2 kpc (adjusted to a distance of 785 kpc), identical to Model N. Seigar et al. (2008) obtained R_d = 5.91 ± 0.27 for the IRAC 3.6 μm profile, which compares well with our Model F.

In the spirit of McD09, we also explored the effect of spiral arm clipping to recover the underlying old stellar population disk, as discussed in Section 4. Such clippings typically yield a shallower scale length, at least for 1D SB profiles (Section 4).

Given the very large number of disk pixels and negligible bulge light into the disk beyond R_min ≃ 1.2 kpc, the disk ellipticity is well constrained by our different analyses (see Tables 2 and 3). Whether through isophotal fits, MCMC or GALFIT, we find for the IRAC data _d = 0.73 ± 0.01. This ellipticity corresponds to a dust-free, razor-thin disk projected at a 74°angle. Assuming a thickness for the M31 disk corresponding to c/a = 0.175 (Haynes & Giovanelli 1984), the M31 viewing angle is then 776.

As we shall see in Figure 9, the disk extends out to roughly R_min ∼ 15 kpc or R_maj ≃ 56 kpc. Other lines of evidence (e.g., Ibata et al. 2005; Worthey et al. 2005) point to the M31 stellar disk extending out to R_maj ∼ 50 kpc, with occasional detections of features (not spatially continuous) with disk-like kinematics out to R_maj ∼ 60 kpc (Ibata et al. 2005). The latter is far beyond the disk–halo transition at R_min ∼ 8 kpc (see Figure 9).

Zoom In Zoom Out Reset image size

The deep star counts accentuate the presence of the halo at R_min > 10 kpc that is seemingly neither the extension of the inner low-n bulge or the exponential disk, nor the manifestation of a faint de Vaucouleurs halo (Tanaka et al. 2010 fit their extended Subaru Suprime-Cam surface brightnesses of the M31 stellar halo with a Hernquist profile and scale radius equal to 14 kpc). As we saw in Section 6, the Sérsic Model W was rejected on account of kinematic constraints. Model U, shown in Figure 8, appears to be a better representation of the halo star counts. Note that the disk and faint components are equal at R_min ∼ 8–9 kpc.

The star counts are consistent with a halo structure best described with a 2D spatial density distribution power-law index −2.5 ± 0.2 (−2 × α as defined in Equation (6)). In three-dimensional, the halo can thus be approximated as a sphere with log ρ(r) = − 3.5log r. The solution from Model U agrees well with other power-law models of the M31 stellar halo (Irwin05; Guha05; Ibata et al. 2007; G09; Tanaka et al. 2010). Tanaka et al. (2010) obtained deep Subaru imaging of M31 along both minor axes, north and south, to find a 2D power-law index equal to −2.17 ± 0.15, in excellent agreement with our value. Our combined results are marginally consistent with Ibata et al. (2007), who measured −1.91 ± 0.11 for their halo power-law fit.

A comparison with the Milky Way is worthwhile. For instance, Keller et al. (2008) modeled the number density of RR Lyrae stars along the ecliptic over the range 10 < R < 45 kpc as a power law with index n = −2.48 ± 0.09 (see also Miceli et al. 2008). Carollo et al. (2010) also identified inner- and outer halo components of the Galaxy based on star counts, kinematics, and abundances of 32360 SDSS and SDSS-II stars. The transition between each component would be around 15 kpc (along the minor axis). The slopes of the power-law spatial density distributions of their inner and outer halo are −3.17 ± 0.20 and −1.79 ± 0.29, respectively. The density profile for their measured outer halo is substantially shallower than that of their inner halo on account of the higher velocity dispersion and higher flattening of the former. Overall, the Milky Way outer halo stars possess a wide range of orbital eccentricities, exhibit a clear retrograde net rotation, and are drawn from a metallicity distribution function that peaks at [Fe/H] = −2.2 (Carollo et al. 2010).

Considering its transition range and power-law slope, our halo component for M31 seems analogous to that of the Milky Way outer halo component identified by Carollo, Kenner and others. Gilbert et al. (2007, 2009b) find the velocity dispersion of M31 halo stars in the range 30–60 kpc in excess of 100 km s⁻¹. This is not only consistent with halo kinematics, but also matches the values quoted by Carollo et al. (2010) for the Milky Way outer halo.

Comparisons between the Milky Way and M31 are still challenging considering the few available lines of sight for the measurements of M31 kinematics (Gilbert et al. 2007, G09). As a result, testing for the rotation of M31's halo is still beyond reach. However, the M31 halo appears to be more metal-rich (〈[Fe/H]〉 = −0.47 ± 0.03) closer to its center (around R_min = 11–20 kpc), and more metal-poor farther out (〈Fe/H〉 = −0.94 ± 0.06 at R_min ∼ 30 kpc, and 〈Fe/H〉 = −1.26 ± 0.1 at R_min > 60 kpc) (e.g., Kalirai et al. 2006; Koch et al. 2008). This qualitatively resembles the findings from Carollo et al. (2010) that the Milky Way's "inner halo" is more metal-rich than its "outer halo." Their "inner halo" (R_min = 10–15 kpc) has a metallicity distribution that peaks at 〈Fe/H〉 = −1.6, while that of their "outer halo" (R_min = 15–20 kpc) peaks at 〈Fe/H〉 = −2.2. Note that the stars in their Milky Way "outer halo" component are within the radial range of the innermost M31 fields of Gilbert et al. (2007, G09). The metallicity of M31's halo at large radii is still quite uncertain: whether it becomes increasingly metal-poor with radius or continually decreases out to, say, R_min ∼ 50 kpc and then levels off is still a point of contention.

Numerous studies of the Milky Way halo based on the distribution and kinematics of RR Lyrae, blue horizontal branch, and giant halo stars over different radial ranges report a fraction of stellar halo light to total galaxy light of ∼1% (e.g., Morrison 1993; Chiba & Beers 2000; Purcell et al. 2007; Bell et al. 2008). This number grows to 2% if unbound Sagittarius stream stars are included (Law et al. 2005). Considerable Milky Way halo light thus comes from dynamically young material in various tidal streams. By contrast, our power-law models for the M31 halo (Table 4, Models R–U, and Equation (7)) yield a fraction of halo-to-total light of ≲ 4%. This value is an upper bound as we assume the continuation of the halo power-law density well into the bulge and disk of M31. The total light of the halo is integrated out to 200 kpc.

Purcell et al. (2008) studied the accretion histories of systems ranging from small spiral galaxies to rich galaxy clusters. Figure 10 shows the predicted light fraction of the diffuse halo stars against the total galaxy light for a variety of accretion histories which account for much of the inner halo light, as a function of host halo mass. Overplotted with blue and red arrows are the measured light fractions, ∼2% and ∼4%, for the Milky Way and M31, respectively. The masses for the Milky Way (≲ 10¹² M_☉) and M31 (∼1.3 × 10¹² M_☉) are taken from Courteau & van den Bergh (1999).¹⁴ Our results are thus consistent with a common evolution of the old stellar systems in the Galaxy and the Milky Way: the inner regions (R < 20–25 kpc) would contain both accreted and in situ stellar populations while the outer regions (R > 25–30 kpc) would be assembled through pure accretion and the disruption of satellites (Carollo et al. 2007, 2010; Purcell et al. 2008 and references therein).

Zoom In Zoom Out Reset image size

Using Equations (3), (5), (7), and (8), we can compute the light fractions of the Sérsic bulge, exponential disk, and power-law halo at any radius as shown in Figure 9 for Models U and W. We also first compute the nucleus total light fraction using KB99's exponential model (Section 3.1). We find that the nucleus contributes a mere 0.05% of the total galaxy light out to 22 kpc (extent of the IRAC 3.6 μm 1D profile along the major axis). We will ignore the nucleus component, and its negligible contribution to the total luminosity profile, in measurements below.

Using the IRAC 3.6 μm 1D brightness profiles, we find that the B/D-only fit yields ∼21% of the bulge light and 79% of the disk light, for a B/D ratio of 0.27 out to R ≃ 22 kpc. Our 2D decompositions of the IRAC image (see Section 4) yield light fractions for the bulge and disk equal to 29% and 71%, respectively, for a B/D ratio of 0.41.

A halo fraction can be estimated from the decomposition of the composite profiles; see Tables 4 and 5. Using Model U, we find that the bulge, disk, and halo each contribute roughly 23%, 73%, and 4% of the total light of the galaxy out to 200 kpc along the minor axis. These light fractions are roughly similar to those obtained with B/D-only fits within 23 kpc along the major radius; this is because the amount of additional B/D light beyond that radius is small and the halo profile takes over at that point.

Figures 8 and 9 show that the light fraction of the galaxy components are consistent with a dominant bulge within R_min ≲ 1 kpc, a dominant disk between 1.5 ≲ R_min ≲ 8 kpc, and a dominant halo beyond R_min ≳ 9 kpc. Note that the latter result is marginally consistent with the light fractions derived by Kalirai et al. (2006), Gilbert et al. (2007), and G09 based on their measurement of stellar kinematics along the minor axis of the M31 halo. Their results show that a disk is undetected as a cold component on the minor axis at R_min ∼ 9 kpc. Taken at face value, the Gilbert et al. constraint implies that much of the light at 9 kpc ascribed to the disk in our decompositions must come from either the bulge or faint components.

Subtle differences between our photometric and spectroscopic analyses could be reconciled if various factors are accounted for. For one, our photometric and spectroscopic data were extracted in different ways and may weigh the stellar density field differently. For instance, our light profile wedges do not sample the same lines of sight as those for the spectroscopically mapped star counts.

The photometrically (this paper) and kinematically based decompositions may also be sensitive to the M31 structure in different ways. The disk structure at R_min ≃ 9 kpc (roughly 35 kpc along the major axis) in the disk plane may be relevant. The disk inclination and position angle are assumed to be constant to these large radii and the kinematic analysis assumes a cold disk. However, at 5–6 disk scale lengths, disks can warp, flare, etc. Resolution of any discrepancies between our respective decompositions ultimately relies on a simultaneous analysis of the photometric and spectroscopic data.

Our paper has provided an investigation of luminosity profile decomposition methods. Whether surface brightness data are (1) binned or not, (2) treated in magnitudes or counts, (3) clipped for non-axisymmetric features, (4) analyzed in 1D or 2D, (5) originate from wedges (cuts) or azimuthally averaged profiles, and/or (6) have different radial ranges or are modeled using different minimization schemes (e.g., MCMC versus NLLS), they may yield slightly different structural parameters. The uniqueness of model parameters in bulge+disk+halo decompositions is also often thwarted by parameter covariances. Despite these caveats, we have achieved self-consistent decomposition results. The technical considerations about light decomposition methods presented in this paper (see also Mac03) suggest systematic errors of galaxy structural parameters from method to method (and author to author) of the order of 20%.

Perhaps the most significant difference between 1D and 2D model parameters (fitted over the same radial range) comes from the ability to decompose the position angles and ellipticities of two photometrically distinct superimposed systems. This crucial operation is only possible in 2D. Our simultaneous fitting of the minor and major axes in 1D has also provided independent information about B/D ellipticities, but not position angle. Future investigations of the M31 2D light distribution will require extensive MCMC decompositions, which are currently computationally expensive (e.g., Yoon et al. 2011).

The extended M31 galaxy light profile with deep star counts reveals a halo component beyond R_min ≃ 10 kpc at a brightness level μ_I ∼ 27 mag arcsec⁻² (Figure 5; also Guha05; Ibata et al. 2007; Tanaka et al. 2010). Because new, deep galaxy surveys¹⁵ should reach at least a few orders of magnitude below this brightness threshold, the common detection of galaxy halos should soon be within reach. If M31 is at all representative, this is important for galaxy structure studies since the decomposition of deep galaxy light profiles will likely be skewed if the halo component is ignored.

M31's n = 2.2 bulge is mostly dominant over the range R_min < 1.2 kpc. The disk takes over in the range 1.2 kpc <R_min < 9 kpc, whereas the halo likely dominates the light for R_min > 9 kpc (see Figure 8). Our intimation of disk dominance between 1.2 and 9 kpc on the minor axis is only marginally consistent with stellar kinematics of the M31 halo whereby the disk is completely cold and absent at 9 kpc (Gilbert et al. 2007; G09). We stress that our decompositions are based on photometric information alone; additional information about the overall structure of M31, such as radial abundance variations or the complete kinematic mapping of M31's major structural components, may possibly yield somewhat different decompositions than the ones presented here. Any tension that may exist between the structural decompositions derived from our photometric analysis and separate kinematic analyses can only be resolved with a simultaneous modeling of the two data sets. Such an exercise is beyond the scope of the present paper but will undoubtedly drive future investigations of the M31 galaxy structure.

The ratio of bulge-to-disk scale lengths, R_e/R_d, has been viewed as a telltale signature of galaxy's secular evolution (Courteau et al. 1996; Mac03; Kormendy & Kennicutt 2004). Mac03 showed that most nearby spiral galaxies have a B/D ratio R_e/R_d ≃ 0.2, which is a natural prediction of secular evolution models (Courteau et al. 1996). It is thus remarkable that, for M31, R_e/R_d ≃ 0.2. Furthermore, a value of n < 2, as for M31, is often interpreted as evidence for secular evolution (Courteau et al. 1996; Fisher & Drory 2010). M31's R_e/R_d ratio and bulge Sérsic index close to two suggests some amount of secular evolution for the bulge.

However, a discussion of the formation scenarios of the M31 bulge cannot simply rely on the shape of its luminosity profile. For instance, the latter may be biased by a recent frosting of star formation (e.g., de Jong & Lacey 2000; MacArthur et al. 2004). The detailed investigation of bulge formation scenarios requires resolved spectroscopy. MacArthur et al. (2009) used deep Gemini/GMOS spectra to show that numerous late-type spirals with nearly exponential (n ∼ 1) bulges harbor mostly old stellar populations. This is confirmed for M31 by Saglia et al. (2010), whose long-slit spectra centered on the M31 bulge show a predominantly old (⩾12 Gyr), slightly α-element enhanced ([α/Fe] ≈ +0.2) solar metallicity bulge. In the inner arcseconds, the M31 luminosity-weighted age drops to 4–8 Gyr and the metallicity increases to three times that of solar. A possible formation scenario thus involves a very rapid, early assembly of the M31 bulge followed by subsequent accretion via secular transport of enriched, disk material to the center. This is neither a purely classical nor secular evolution formation scenario, but rather one where both mechanisms are at play. Similar conclusions have been drawn for the Milky Way; Babusiaux et al. (2010) report the presence of two distinct populations in Baade's Window with a metal-rich population with bar-like kinematics and a metal-poor population with kinematics indicative of an old population. The entangled interpretation of spiral galaxy bulges, like M31's, in terms of both merger-induced (classical) and secular evolution for larger samples of galaxies is also discussed in Peletier et al. (2007), MacArthur et al. (2009), and Roediger et al. (2011), among others.

Numerous questions remain about the spatial distribution, populations, and kinematics of M31 stars. However, one conclusion is clear: M31 likely evolved through a rich past that involved early "classical" merging and structural readjustments via secular evolution effects and satellite accretion. As a result, M31's luminosity profile does not follow a simple de Vaucouleurs law.