The Complex Structure of the Bulge of M31

The Sérsic function is defined as

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{e}\exp \left[-\kappa \left({\left(\displaystyle \frac{r}{{r}_{e}}\right)}^{1/n}-1\right)\right].\end{eqnarray} \tag{ 1 }$$

The Sérsic function has seven parameters, including the elliptical parameters: x_c , y_c , m_tot, r_e , n, q, θ_PA, with m_tot the integrated magnitude, r_e the effective radius, and n the Sérsic index. The exponential profile has the form

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}\exp \left(\displaystyle \frac{-r}{{r}_{s}}\right),\end{eqnarray} \tag{ 2 }$$

with six parameters: x_c , y_c , m_tot, r_s , q, θ_PA. The Gaussian profile has the form

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}\exp \left(\displaystyle \frac{-{r}^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 3 }$$

with six parameters: x_c , y_c , m_tot, FWHM, q, θ_PA, with FWHM = 2.354σ. The Moffat profile is given by

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)=\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{{\left[1+{\left(r/{r}_{d}\right)}^{2}\right]}^{n}},\end{eqnarray} \tag{ 4 }$$

with ${r}_{d}=\tfrac{\mathrm{FWHM}}{2\sqrt{{2}^{1/n}-1}}$. The seven parameters are x_c , y_c , m_tot, FWHM, concentration index n, q, and θ_PA. The Nuker profile takes the form

$$\begin{eqnarray}&&I(r)={I}_{b}{2}^{\tfrac{\beta -\gamma }{\alpha }}{\left(\displaystyle \frac{r}{{r}_{b}}\right)}^{-\gamma }{\left[1+{\left(\displaystyle \frac{r}{{r}_{b}}\right)}^{\alpha }\right]}^{\tfrac{\gamma -\beta }{\alpha }}.\end{eqnarray} \tag{ 5 }$$

The Nuker profile has a total of nine free parameters: x_c , y_c , μ_b , r_b , α, β, γ, q, θ_PA, with μ_b the surface brightness at radius r_b . Additional available functions are described in Peng et al. (2002) and the GALFIT documentation on the GALFIT web page. ³

To account for telescope optics and, for ground-based telescopes, atmospheric seeing, GALFIT convolves the model image with a PSF image. To extract the PSF image for each UVIT filter, the following process was implemented. For the UVIT image in each filter, several isolated point-source stars were identified. From these we chose one where the star was best centered in the 5 × 5 array of pixels around the star. CCDLAB produces images with 04168 pixels, which was designed to oversample the instrument resolution of ≃1'' by a factor of ∼2.5. Because the PSF has extended wings (Tandon et al. 2017), we selected a subimage of 55 × 55 pixels centered on the point source. Because there is diffuse light in M31, e.g., from the bulge and disk, we subtracted the mean of the four 7 × 7 pixel corners from the image to produce the diffuse background-subtracted PSF. The GALFIT modeling for each UVIT filter image of M31 utilized the corresponding filter's PSF fits image.

To model as much of the bulge light as possible while avoiding the light from the disk and inner spiral arms, the central 1001 × 1001 pixel (4172 × 4172) area of Field 1 was chosen. We visually inspected larger regions (1101 and 1201 pixels) but could see some light from the spiral arms in those images at the corners, so we rejected use of larger regions. The image of the chosen region is shown in Figure 1 for the five UVIT FUV and NUV filters: F148W, F169M, F172M, N219M, and N279N.

Zoom In Zoom Out Reset image size

An initial array of GALFIT tests was performed with the goal of determining the approximate fit qualities of the different radial profile functions of GALFIT (we tested the five functions listed in Section 2.2.1) for both the bulge and nucleus using the F148W image. The model included the three brightest stars in the chosen 1001 × 1001 pixel area. ⁴ The Sérsic, Gaussian, Moffat, Nuker, and PSF functions in GALFIT were used as model components in trials where up to two components were used to fit the bulge region. All model components were elliptical (Peng et al. 2002) because galactic components (e.g., bulge, disk, and halo) projected onto a two-dimensional image are elliptical in shape.

Next, we carried out model fits to test the constancy of the model parameters versus wavelength. Three separate series of tests with a Sérsic bulge and a Moffat nucleus were performed with (i) all free parameters, (ii) a constant effective radius, and (iii) a constant Sérsic index. After performing an analysis of the parameters versus wavelength for each of the three tests, we conclude that none of the parameters are consistent with being constant versus wavelength. Thus, further modeling requires Sérsic index, FWHM, or concentration index as free parameters for each UVIT filter.

Dust cloud absorption features are visible northwest of the center of the bulge in the UVIT images (see the panels for the various filters in Figures 1(a)−(e)). We excluded the areas with strong dust cloud features from the fits by using a mask because our purpose was to fit the bulge structure and not the dust clouds. We used a Python script to generate a list of pixels that correspond to the coordinates of the pixels affected by the main dust cloud regions visible in the bulge. The F148W image was used to identify the dust clouds owing to its better S/N compared to the images taken with the other four filters. The mask is shown in Figure 1(f) and was used for the model fits for all five UVIT filters.

In addition to the stellar light from the bulge, there are several bright foreground stars in the image. We test adding a different number of stars in the model for the bulge. Rows 4−6 of Table 2 show the χ² for model fits including the bulge (Sérsic plus Gaussian) plus one, two, or three stars for comparison to the χ² for a bulge model with a Sérsic plus Gaussian and no stars (row 3). Including the brightest star, the two brightest stars, or the three brightest stars successively improves the χ² by ∼15,000, ∼1500, and ∼4500 compared to no stars. ⁵

We did not include more stars because the improvement in fit quality (χ²) was significantly less than for the three brightest stars. For all of the fits presented here the background and the positions and magnitudes of the three brightest stars are free parameters (three per star), to give the number of background plus star parameters of 10.

The nucleus of M31 has a complex structure as shown by Peng (2002), which includes the compact nucleus and the nuclear bulge surrounding it. The compact nucleus was resolved into double components by Lauer et al. (1993) with Hubble Space Telescope (HST) Wide Field and Planetary Camera (WFPC) observations. The HST WPPC2 F555W (V-band) image was analyzed by Peng (2002) using GALFIT two-dimensional fitting of the central 5'' × 5'' area. Kormendy & Bender (1999) show that the double nuclei in V band (called P1 and P2) are the two ends of an eccentric disk of stars orbiting the central black hole, which has mass ∼3 × 10⁷ M_⊙. Peng (2002) finds that the central region is well modeled by four physical components: P1, P2, the M31 elliptical (q ≃ 0.81) bulge, and a spherical nuclear bulge (R_e ≃ 32).

In order to fit the brightness distribution of the bulge only, we exclude from the fit the emission associated with the nucleus (compact nucleus plus nuclear bulge). For this purpose a series of GALFIT pixel masks of square shape and varying size, centered on the nucleus, were created. The best-fit χ² was given by the nucleus pixel mask of dimensions 41 × 41 pixels (17'' × 17''). Table 3 shows the results of tests conducted to verify whether masking the nuclear region improves the fit of the bulge. The χ² is reduced at high significance (several tens of σ; Press et al. 1992), ⁶ indicating that there is a separate component (the nucleus) smaller than the 41 × 41 pixel mask. The fact that the bulge region is fit significantly better with the nuclear mask for the tested cases (one-Sérsic, two-Sérsic, and three-Sérsic functions) is a strong indicator that the nuclear region has additional components, in agreement with the results of Peng (2002). Hereafter, when fitting the "bulge only," the nucleus mask is included to omit the nuclear region (see Section 3.1). When fitting the bulge plus nucleus, the nucleus mask is omitted, to include the nuclear region (see Section 3.2).

For these fits, the Sérsic function for the nuclear bulge has effective radius R_e and index n fixed to the values from Peng (2002). We use a PSF function for the compact nucleus because it is unresolved with UVIT, although it is complex at the subarcsecond scale (Peng 2002). Table 6 shows results (rows 5−12) for model fits including the nuclear region (no nucleus mask). These use various numbers of Sérsic functions for the bulge, plus one Sérsic for the nuclear bulge and a PSF function for the compact nucleus.

Comparing the different models in Table 6, we find the following. First, we compare fits for the bulge only (rows 1−3 in Table 6) with the fits including nucleus (row 5) or nuclear bulge (rows 6−8), with the nuclear bulge parameters fixed to the (Peng 2002) values, except normalization. For fits including the pixel mask, the degrees of freedom (dof) is reduced by 1681, from 1,002,001 − N_par to 1,000,320 − N_par. The three-Sérsic plus PSF fit (row 5) for the bulge plus nuclear region (with 1681 more pixels than the bulge with nuclear mask) has higher χ² by ∼6200 than the three-Sérsic fit for the bulge alone (row 3). This is significantly higher than the expected increase of 1681, the extra dof introduced by including the 41 × 41 pixel nuclear region. This shows that the PSF alone does not provide a good fit to the nuclear region. Similarly, comparison of row 1 with row 6, row 2 with row 7, or 3 with row 8 shows increases of χ² of ≃5000 for each case. This is much higher than the extra dof introduced, showing that the nuclear bulge with fixed parameters does not provide a better fit.

Next, we compare fits to the bulge only (rows 1−4) with fits including bulge, nucleus, and nuclear bulge (rows 9−12). The one, two, three, and four Sérsic bulge fits with nucleus and nuclear bulge have increases in χ² of 3841, 3319, 4167, and 3724, respectively, over fits without nucleus and nuclear bulge. This is much greater than the extra dof introduced. We conclude that the model with fixed parameter Sérsic plus PSF is not a good model for the nuclear region (nucleus plus nuclear bulge).

For our final set of fits to the bulge plus nuclear region, we use models for the compact nucleus plus nuclear bulge, with nuclear model parameters not tied to the values from Peng (2002). The nuclear region observed by UVIT in the F148W filter has lower resolution than the HST F555W data, and the wavelength is significantly different, so different stars are prominent in F148W than in the HST data. The F148W filter extends from 125 to 175 nm (Tandon et al. 2017), so the bulge and nuclear region structure can be quite different than what was found at 555 nm by Peng (2002).

The procedure for finding the best-fit model, without fitting beyond what is statistically significant, is outlined in Figure 2. We start with a single function to model the bulge, including the nuclear region, and compare it to fits with the other four functions (see Section 2.2.1). The function with the best fit is chosen, and then the complexity is increased by adding an additional function to the function set. All five functions are tested as the added function, and the function set with the best fit was kept. Next, we test whether the current function set has a significantly better fit than the best fit with the function set with one fewer function. If yes, then we test a function set with one additional function, as shown in the flowchart. If not, the added function is not justified by the data, and the set with one fewer function is the final best-fit function set.

Zoom In Zoom Out Reset image size

Models with free bulge and nuclear parameters with four to eight components are shown in Table 7, with parameter values and χ² values given. The model with nine components does not yield a statistically better fit and so is not shown. Neither do we show the models with three or fewer components, as these are inadequate to fit the bulge plus nuclear region and have significantly worse fits. The parameter errors are those generated by GALFIT from the fitting process.

We compare the models to comparable models with fixed nuclear parameters, shown in Table 6. Model 1 in Table 7 has a lower χ² by 2479 than Model 8 from Table 6 but only two more parameters, thus a significantly better fit. Comparing Model 1 in Table 7 with Model 11 from Table 6 yields that Model 1 in Table 7 has a lower χ² by 1268 with one less parameter. Thus, we conclude that the free nuclear region parameters result in much better fits and that the nucleus shape for the 148 nm image is significantly different than for 555 nm.

We next compare the different models in Table 7. The models with more components (and thus parameters) generally have lower χ². The Δχ² values are compared here to determine how much better are the more complex models. Model 2, with 41 parameters, has Δχ² = 815 compared to model 1, with six extra parameters. This is a highly significant improvement in fit (Press et al. 1992). Model 3 has Δχ² = 1351 compared to model 2, with seven extra parameters, again a highly significant improvement in the fit. Model 4 has Δχ² = 123 compared to model 3, with six extra parameters. This is a significant improvement in the fit. ⁸ Model 5 has Δχ² = 69 compared to model 4, with five extra parameters, which is a significant improvement in fit. We did not proceed with more complex models than Model 5 for two reasons: the improvement in fit was no longer highly significant, and GALFIT would not converge from almost all the models attempted, likely because the number of parameters was too large. Thus, our final model for the bulge and nuclear region is Model 5 with three components for the bulge and five components for the nuclear region.

Residual images are obtained by subtracting model images from the original image. Figure 3 illustrates the process of finding the best combination of functions for the bulge of M31. It shows the original F148W observed image (top left) and, in the remaining panels, residual images for successively better models for the central 125'' × 125'' (475 pc × 475 pc) region centered on the M31 compact nucleus.

Zoom In Zoom Out Reset image size

The best-fit χ² values are successively lower for the residual images in successive panels. The second panel in the top row is for a single-Sérsic model for the bulge and no nuclear model (row 1 in Table 6). The single-Sérsic model, although it accounts for the bulk of the light from the bulge (in comparison with the first panel), has a significant residual in the nuclear region (inner 41 × 41 pixel region) and thus does not fit the light distribution of bulge plus nuclear region. Successive models (residuals in panels 3 and 4, for the two-Sérsic and three-Sérsic models) fit the region outside the 41 × 41 pixel nuclear region and represent the structure of the outer region more accurately, but they do not fit the nuclear region. Panels 5−9 are for models 1−5 in Table 7, which include components for the nuclear region. These show increasingly better fits (smaller residuals) in the nuclear region while maintaining or improving the fit to the outer region. The last panel in the bottom row is for our final bulge plus nuclear region model 5 with eight components shown in Table 7. This model yields no visible artifacts in the nuclear region at the center or in the bulge farther out.

From the residual images in Figure 3, it is seen that as the fit improves, the residuals in the nuclear region decrease. The final model fits well the large-scale structure of the bulge and the smaller-scale structure of the nuclear bulge and nucleus. There remains significant small-scale structure in the last image shown. The residual structures that are not modeled are seen to consist primarily of dust lanes and a number of small (few arcseconds across) bright spots scattered throughout the region.

The two bulge categories of pseudobulges and of classical bulges are essential in understanding galaxy formation history (KK04; Fisher & Drory 2008). According to KK04, the Sérsic index can be used to classify bulges, or at least provide evidence in favor of one of the types in the absence of kinematic information. A Sérsic index value less than 2 is consistent with a pseudobulge, while classical bulges have a characteristic Sérsic index greater than 2 (KK04; Gadotti 2008; Fisher & Drory 2008, 2016). Athanassoula et al. (2016) distinguish (i) classical bulges by having a spheroidal shape, Sérsic index >2.5, and ${V}_{\max }/\sigma$ values that are less than for isotropic oblate rotators; and (ii) pseudobulges by having the shape of a disk, Sérsic index <2.5, and clear signs of rotation.

Caon et al. (1993) demonstrated that ellipticals and galactic bulges are better fit if the exponent of the r^1/4 term in the de Vaucouleurs formula is a free parameter as r^1/n , equivalent to the Sérsic function with index n. Further studies, such as Courteau et al. (1996), strengthened the idea that "late-type" bulges have a Sérsic index between 1 and 2. This led KK04 to include n between ≈1 and 2 as one of the diagnostic tests for galactic pseudobulges. Fisher & Drory (2008) performed bulge–disk decompositions on a sample of spiral galaxies and likewise concluded that Sérsic indices smaller than 2 are indicative of a pseudobulge, while those greater than 2 are indicative of a classical bulge. This is confirmed by the Gadotti (2008) study on the structural parameters of bulges, bars, and disks. They found that bulges show significantly different parameters than ellipticals, with a greater dispersion for pseudobulges. Fisher & Drory (2010) find that values of n < 2 may be interpreted as evidence of secular evolution in the bulge, and they suggest that their n value close to 2 suggests "some amount" of secular evolution. Nine S0–Sb galaxies that have two-component bulges (classical plus pseudobulge) are the subject of an observational study by Erwin et al. (2015). They find that the disky pseudobulges consist of an exponential disk with nuclear rings, nuclear bars, or spiral arms and that the classical bulge components have Sérsic indices in the range 0.9–2.2.

Simulations of galaxy evolution provide a guide to the structural differences between pseudobulges and classical bulges. Guedes et al. (2011) present a cosmological N-body hydrodynamic simulation in which the evolutionary pathway of a close analog to the Milky Way disk galaxy is studied. Bulge formation was observed to have occurred through slow secular processes by the interaction of the various structural components of the galaxy, rather than fast and energetic merger events, which were absent in the simulation. Their simulation outcome of a "present-day" disk shape with a complex structure showing nuclear, bulge, disk, and halo regions was subsequently subjected to a GALFIT decomposition. The resulting bulge Sérsic index was n = 1.4. Simulations for a major merger of two gas-rich pure disk galaxies (Keselman & Nusser 2012) find that the resulting bulge has Sérsic index ∼0.3–0.9. When the initial galaxies have gas halos (Athanassoula et al. 2016), the immediate outcome is a classical bulge, but a pseudobulge and a thin disk grow from gradually accreted gas, which can be considered secular evolution. The kinematics and range of stellar ages for the classical and pseudobulge are very different. The classical bulge (stars born at time t < 1.4 Gyr) has a Sérsic index of n ∼ 4–6; the four remaining components (stars born at time intervals 1.4 Gyr < t < 1.8 Gyr, 1.8 Gyr < t < 2.2 Gyr, 2.2 Gyr < t < 9 Gyr, and 9 Gyr < t < 10 Gyr) can be fit with a combination of three exponential disks.

A second characteristic of pseudobulges as described by KK04, Courteau et al. (1996), and Fisher & Drory (2010) is the presence of boxiness. Bulges showing box-shaped isophotes are well established observationally and theoretically (Sandage 1961; de Vaucouleurs 1974). Boxy bulges are observed in about one-fifth of edge-on spiral galaxies (de Souza & Dos Anjos 1987; Lütticke et al. 2000; Laurikainen & Salo 2016). KK04 conclude that the detection of boxy bulge isophotes is sufficient for the identification of a pseudobulge. This stems from the connection boxy bulges seem to have with galactic bars: they appear to be created by the interaction of the bar with the disk, as secular evolutionary mechanisms slowly construct the boxy bulge out of disk material.

For M31, Courteau et al. (2011) use one-dimensional decomposition methods to model the bulge and obtain a Sérsic index n = 2.2 ± 3. They measured a ratio of bulge scale length to disk scale length of 0.2, which agrees with the value of 0.2 that is obtained from secular evolution models (Courteau et al. 1996). Courteau et al. (2011) note that their Sérsic index is close enough to 2 to be consistent with some amount of secular evolution in the bulge, according to Courteau et al. (1996), Fisher & Drory (2010), and KK04. The study here shows that a single-Sérsic fit to the five UVIT images (Table 8) yields indices between ∼2.1 and ∼2.6. Here we used two-dimensional image fitting with GALFIT, which should better model the light distribution than one-dimensional fits (see Section 4.1). We obtain an M31 bulge Sérsic index of 2.14 for the 279 nm UVIT image, consistent with earlier results, and indices of ≃2.5 for the four shorter-wavelength UVIT images (148–219 nm). The fact that the Sérsic index is ≃2–2.5 and neither clearly <2 nor clearly >2.5 indicates that there are probably contributions from a classical bulge and a pseudobulge in M31. The four shorter-wavelength images yield a larger Sérsic index, possibly indicating a smaller contribution from a pseudobulge relative to the classical bulge at shorter wavelengths.

Several studies indicate that the M31 bulge is intrinsically an oblate spheroid with an axis ratio of ∼0.8 (Kent 1983; Widrow et al. 2003; Geehan et al. 2006), while others indicate a triaxial nature (Lindblad 1956; Stark 1977; Kent 1989; Stark & Binney 1994; Berman 2001; Berman & Loinard 2002). For both cases, two-dimensional fitting of the image is superior to one-dimensional fitting.

A near-infrared mapping study of M31 conducted by Beaton et al. (2007) revealed the boxy morphology of the bulge with high S/N and with little interference from dust extinction. The UVIT FUV and NUV images of the bulge of M31 show a boxy morphology. Figure 5 shows contour lines from the five UVIT images (panels (a)–(e)). Panel (f) shows the contour lines from Beaton et al. (2007), which illustrate similar boxiness to that seen in the N279N image in panel (e). The other UVIT filter images exhibit boxiness, although less strongly than that seen for the N279N image.

Zoom In Zoom Out Reset image size

Most of the boxy or peanut-shaped features that have been reported in galaxies have been detected in infrared bands (Kormendy & Barentine 2010; Laurikainen & Salo 2016), such as Spitzer 3.6 μm images. In the UV there have been no previous reports of bulge boxiness. However, all five of the UVIT filter images of the bulge of M31 show signs of a boxiness. Considering KK04's third criterion, the bulge boxiness indicates a pseudobulge interpretation for M31's bulge.

There are two ways to measure diskiness or boxiness from fitting images using GALFIT. The first method is the use of generalized ellipses instead of ellipses, as follows. The radial functions are elliptical in shape, accomplished by making the radial coordinate r constant on ellipses, defined by

$$\begin{eqnarray}&&r(x,y)={\left(| x-{x}_{c}{| }^{2}+| \displaystyle \frac{y-{y}_{c}}{q}{| }^{2}\right)}^{\tfrac{1}{2}}.\end{eqnarray} \tag{ 6 }$$

To allow the ellipse to have arbitrary orientation, (x, y) are rotated coordinates. To change this into generalized boxy or peanut-shaped ellipses, the power is changed from 2 to C_o + 2 as follows (Athanassoula et al. 2016):

$$\begin{eqnarray}&&r(x,y)={\left(| x-{x}_{c}{| }^{{C}_{o}+2}+| \displaystyle \frac{y-{y}_{c}}{q}{| }^{{C}_{o}+2}\right)}^{\tfrac{1}{{C}_{o}+2}}.\end{eqnarray} \tag{ 7 }$$

C_o is the parameter that controls the diskiness or boxiness of the isophotes (lines of constant r(x, y)). For C_o = 0 the isophotes are perfect ellipses. For C_o < 0 the isophotes becomes more disk-like or disky as C_o becomes more negative. For C_o > 0 the isophotes become more boxy as C_o becomes more positive. Table 9 shows results of fitting generalized ellipses for the five UVIT images. All values of C_o are >0, showing that the M31 bulge has a boxy shape.

The second method is to fit radial functions modified by azimuthal functions in GALFIT, as described in Peng et al. (2010). The m = 4 Fourier mode describes whether the bulge is boxy or disky. Table 9 shows the results of fitting the bulge including the m = 4 Fourier mode. For all five UVIT images the boxiness of the bulge is significant (amplitude >0 at 7σ–25σ). The Fourier mode and C_o tests of boxiness are consistent for each image: boxiness is strongest for the F148W image, less strong for the F172M image, even less strong for the F169M and N279N images, and weakest but still significant for the N219M image.