ASTROSAT/UVIT Near and Far Ultraviolet Properties of the M31 Bulge

AstroSat has surveyed M31 with the UVIT telescope during 2017 to 2019. The central bulge of M31 was observed in 2750-2850 A, 2000-2400 A, 1600-1850 A, 1450-1750 A, and 1200-1800 A filters. A radial profile analysis, averaged along elliptical contours which approximate the bulge shape, was carried out in each filter. The profiles follow a Sersic function with an excess for the inner 8"in all filters, or can be fit with two Sersic functions (including the excess). The ultraviolet colours of the bulge are found to change systematically with radius, with the center of the bulge bluer (hotter). We fit the UVIT spectral energy distributions (SEDs) for the whole bulge and for 10 elliptical annuli with single stellar population (SSP) models. A combination of two SSPs fits the UVIT SEDs much better than one SSP, and three SSPs fits the data best. The properties of the three SSPs are age, metallicity (Z) and mass of each SSP. The best fit model is a dominant old, metal-poor (10^10 yr, log(Z/Zsun)=-2, with Zsun the solar metallicity) population plus a 15% contribution from an intermediate (10^9.5 yr, log(Z/Zsun)=-2) population plus a small contribution (2%) from a young high-metallicity (10^8.5 yr, log(Z/Zsun)=-0.5) population. The results are consistent with previous studies of M31 in optical: both reveal an active merger history for M31.


INTRODUCTION
The spiral galaxy Andromeda, also known as M31, is the nearest such galaxy to our own. Our external vantage point makes it feasible to study aspects of M31 that are difficult to study in our own galaxy. Interstellar extinction is not as major a factor when studying M31 as it is studying the Milky Way, owing to the former's inclination that allows observation of the whole disk of the galaxy. Studies of large numbers of stars have been done in great detail for the Milky Way, but distances often possess high uncertainty, due in part to extinction. An advantage of studying objects in M31 is that it is at a well known distance (785±25 kpc, McConnachie et al. 2005). The intrinsic brightness of many objects can therefore be more precisely measured than galactic sources. M31 has been observed in optical on numerous occasions. The highest resolution observations were carried out with the Hubble Space Telescope (HST), including the Pan-chromatic Hubble Andromeda Treasury (PHAT) survey (Williams et al. 2014). In near and far ultraviolet (NUV and FUV), the GALEX instrument (Martin et al. 2005) has surveyed M31.
The AstroSat mission, launched on September 28th, 2015, has carried out a survey of M31 in near and far UV. AstroSat is an orbiting observatory equipped with five instruments: the UV Imaging Telescope (UVIT) for visible and UltraViolet (UV); the Soft X-ray Telescope (SXT), Large Area Proportional Counters (LAXPC) and Cadmium-Zinc-Telluride Imager (CZTI) instruments for soft through hard X-rays; and the Scanning Sky Monitor (SSM), an X-ray survey instrument (Singh et al. 2014).
Analysis of the M31 UVIT survey observations have been presented in part by Leahy et al. (2018), , ,  and Leahy et al. 2021 (submitted for publication). These papers have dealt with resolved stellar objects, including analysis of UV brightest stars in the bulge, the The images are 410" by 410" across, and centred on 00h 42m 44.34s ra., +41d 16' 08.44" dec. M31 UVIT point source catalog, matching UVIT point sources with Chandra sources in M31, improvements in astrometry and photmetry for the M31 survey, first results from matching UVIT sources with HST/PHAT sources in the NE spiral arms of M31, and a study of FUV variable sources in M31 using a new second epoch observation of the central field of M31. The UVIT observations for the M31 survey are described in the M31 UVIT point source catalog paper . 19 UVIT fields of 28 arcmin diameter were required to cover the M31 area (Fig. 2 of . The current analysis is primarily concerned with the unresolved ultraviolet sources in M31's central bulge. The bulge of a galaxy is the central concentrated set of stars which is ellipsoidal in shape, smoothly distributed and consists of a substantial fraction of the galaxy's stars. UVIT's multiple filters within the NUV and FUV bands allow for an in-depth look at the UV colours of the bulge. In Sections 2 and 3, the observations and data analysis methods are described. The results are reviewed in Section 4 including radial profile analysis and UV colour change of the bulge with radius. In Section 5 synthetic colour magnitude diagrams (CMDs) are used to model the colour change of the bulge with radius and constrain the stellar populations in the bulge. We close with a brief summary.

OBSERVATIONS
The UVIT instrument onboard AstroSat has high spatial resolution ( 1 arcsec), a 28 arcminute field of view, and was capable of observing in a variety of FUV and NUV filters (Tandon et al. 2017). However, partway through observations of M31, the NUV detector of UVIT failed, resulting in partial coverage of M31 in NUV. The entire galaxy has been observed in the FUV, and the central bulge is included in the portion with NUV data. The central bulge of M31 is entirely contained within Field 1 of the UVIT Andromeda survey. The bulge of M31 was observed by UVIT in 2017 (Leahy et al. 2018) and at a second epoch in 2019. New in-orbit calibrations of the UVIT instrument have been carried out by (Tandon et al. 2020) and are applied here to the 2017 and 2019 observations. The basic data processing was carried out using CCD-Lab (Postma and Leahy, 2017), with updated astrometry calibration from , and updated photometry as described in  and Leahy et al. 2021 (submitted for publication). M31 Field 1 (the bulge field) was observed by UVIT instrument in 5 FUV and NUV filters. From longest to shortest wavelength, these are: N279N, N219M, F172M, F169M, and F148W, where the label contains the central wavelength (in nm) of each filter and band 1 .
The images were taken across two epochs, labelled A for the 2017 observation and B for the 2019 observation. The first set of observations included N279N, N219M, F172M, and F148W images, while the second set included F172M, F169M and F148W images. Merged images, here labelled "M", were produced for filters F172M and F148W which had data from both A and B epochs.
1 N is for NUV, F is for FUV. The filter transmission curves are given in Tandon et al. 2017. The merged images were created by adding the counts images for A and B, adding the exposure maps for A and B, then dividing the counts images by the exposure images to create a count/s image from which fluxes can be extracted using the UVIT calibrations (updated in Tandon et al. 2020).
The full 3-colour image of Field 1 made from the 2017 epoch N279N, N219M and F148W images is shown in Leahy et al. (2018). The image includes the bright central bulge, which is likely dominated by unresolved FUV and NUV emission from large numbers of low-luminosity stars in the bulge. Field 1 contains other prominent features, including resolved and unresolved emission from numerous stars in the inner spiral arms of M31, dust lanes which are visible as dark lanes, and several foreground stars, which can be identified using either Gaia parallax  or colour-magnitude diagram analysis .
The central part of Field 1 is the area that contains the central bulge of M31, here defined as the inner 7 arcminute diameter part of the 28 arcminute diameter full area of Field 1. The area studied here is shown in Figure 1, with the images of the bulge shown for the five different NUV and FUV filters.
Due to the low inclination of M31, a few dust lanes are visible, e.g. NW of the bulge by ∼80 arcmin. The amount of area covered by the dust lanes is not large, so we did not correct for the reduction in flux caused by these dust lanes in our analysis. However, we did test an analysis which accounts for the bright point sources in the images. We masked the bright point sources and replaced them with the average values of surrounding pixels. We compared the brightness profiles derived with and without the point sources and found there was no significant difference. Thus we decided to leave the point sources in the images for the final analysis.
3. DATA ANALYSIS 3.1. Elliptical brightness profile extraction Figure 2 illustrates the ellipticity of the central bulge of M31 at the longest NUV and shortest FUV wavelengths. Shape changes of the bulge across different wavelengths are investigated here by fitting elliptical Gaussian functions to the images using CCDLab for a set of square bounding boxes of different sizes all centred on the bulge. The centre of the bulge was chosen as the brightest pixel in the core in the majority of filters. Figure 3 shows the ellipticity (a/b) and inclination (φ) as a function of major axis (a). Inclination is measured in radians clockwise from the EW axis. Data from bounding boxes <40 pixels half-width are not shown. This corresponds to a cutoff of 4.2 , the small central part of the bulge, which has a complex structure (e.g. Lockhart 2017).
To measure the brightness and UV colours of the bulge we measured average flux per pixel as a function of semi-major axis, extracted using elliptical annuli. This method was chosen over radial profiles using circular annuli because the images are not circularly symmetric, and thus the brightness changes significantly with azimuthal angle for circular annuli. The brightness profiles using elliptical annuli accounts for the elliptical shape of the bulge projected on the sky. Pixels (coordinates x,y) were counted within a given ellipse annulus if they satisfied: Here (F 1x, F 1y) and (F 2x, F 2y) are the positions of the focii of the given ellipse: Cy) is the (x, y) position of the centre of the bulge on the image, a is the semi-major axis desired, in pixels, b the corresponding semi-minor axis, and p the inclination of the major axis from the image's horizontal, in radians. The average flux per pixel is calculated by summing the counts/sec in each pixel within a thin elliptical annulus, dividing by the number of pixels in the annulus, then using the UVIT calibration to convert count rate per pixel to flux per pixel.
To extract flux measurements for the same elliptical annuli of the bulge for all filters, we used the ellipse parameters vs. semi-major axis from the N279N A image analysis. As a test, the analysis was repeated using elllipse parameters derived from the F148W M image. However there was little difference in the resulting ellip-tical profiles, so the analysis presented here was carried out using the N279N A ellipse parameters. Average flux per pixel was extracted for 240 elliptical annuli, from the centre out to 400 pixels (167 ) along the major axis. Examples of these ellipses are drawn in figure 2. One ellipse was used for each of the first 200 pixels along the major axis, and then one ellipse for every 5 pixels for the remaining 200 pixels.

UV colour-colour diagram
To study the colour change in the bulge, we produced a colour-colour diagram using the elliptical profiles of flux per pixel from the above analysis. The colours chosen were ratio of flux per pixel of F148W to N219M and ratio of N219M to N279N. The standard deviation of values of the fluxes from the different pixels in each annulus were used as estimate of the flux errors, and propagation of errors was used to calculate the errors in the colours.

Elliptical brightness profile analysis
Analytic functions were fit to the elliptical surface brightness profiles of the bulge to quantify the shape in a simple manner. The fits were carried out using χ 2 minimization. For these fits the fractional error for each Note-a: The uncertainties in all parameters are ∼ 5%.

Note-b:
The fits included data for semi-major axis of 1 to 400 pixels.
Note-c: DOF stands for degrees-of-freedom.
Note-d: The fits included data for semi-major axis of 10 to 400 pixels.
semi-major axis value was taken as 1/ √ N with N the number of counts at tht semi-major axis 2 .
Initially, we fit a Sersic function, f (a), to the observed surface brightness vs. semi-major axis, a. Sersic functions are commonly used to approximate stellar distributions in elliptical galaxies and spiral bulges: with I0 the central surface brightness, n the Sersic index, and k the e-folding constant for decrease of surface brightness. This is equivalent to the original form of the Sersic function given by Sersic (1968), and R e = (b(n)/k) n . For our data, in units of flux per pixel (erg cm −2 s −1 A −1 pixel −1 ), we use units of pixels for a. Figure 5. Double Sersic fits (dashed lines) compared to the data (solid lines). The fits include data from 1 pixels to 400 pixels (167 ). Data errors range from a minimum 1.3 to 1.8% for F148W to a maximum of 2.5 to 4.5% for N279N. Errors for the other filters range between 1.7 to 3.1%. A ±5% error bar is shown on the N219M curve for illustration. The errors are similar to the differences between the data points and the fit curves.
2 The systematic errors in converting counts/s to flux or magnitude (Tandon et al. 2020 andTandon et al. 2017) are small compared to the 1/ (N ) statistical errors for these narrow (1 pixel wide) elliptical areas.
The elliptical profiles generally follow a Sersic shape, with an excess in the center in all filters. The right panel of Figure 4 shows the F148W profile. The Sersic function (black dashed line) fits very well the majority of the bulge (from semi-major axis of 8 to 167 ). The excess is not consistent with a point source at the centre of the bulge, as shown by the red dashed lines in Figure 4, which include Sersic plus the radial profile of a point source (the UVIT point-spread-function or PSF, Tandon et al. 2017). The shape of the excess is too broad and flat compared to the PSF. Higher resolutions are required to study the inner ∼ 2 , as was done by Lockhart (2017) at optical wavelengths. Table 1 shows the results of single Sersic function fits to the bulge, including and omitting the inner 10 pixels (∼ 4 ) of the bulge. Comparing the two cases shows that omitting the inner ∼ 4 results in a significant reduction in χ 2 . Using the fit with the inner 4 omitted, the Sersic index, n, varies with wavelength, with n increasing from 2.3 to 5.7 for decreasing wavelengths 279 nm → 219 nm → 172 nm, then decreasing from 5.7 to 4.1 for 172 nm → 169 nm → 148 nm. The fits are acceptable (χ 2 degrees-of-freedom or DOF) except for the F148W data, which has the highest signal-to-noise of the five filters.
To better fit the F148W data, and to see if the inner part of the bulge could be fit, we carried out fits to the full data (1 to 400 pixels) with double-Sersic functions, and Sersic+Gaussian functions. The Gaussian function is: Results are given in Table 1. For F148W, the Sersic plus Gaussian is significantly better than a single Sersic (χ 2 of 914 vs. 1900), and the double Sersic better than Sersic plus Gaussian (χ 2 of 770 vs. 914). Figure 5 shows the double Sersic fits to measured surface brightness for all filters, illustrating that the central bulge ( 10 ) and rest of the bulge (∼ 10 to 167 ) are well fit. The N219M filter observed surface brightness looks different than the other filters: it is faintest in the inner bulge and approximately matches the shorter wavelength profiles for the outer bulge 3 Figure 6. Left: UV colour-colour diagram of the bulge of M31, for elliptical annuli of different semi-major axes, from 0 to 400 pixels (0 to 167 ). Semi-major axis is encoded by RGB colour, as shown by the panel on the right (the color scale is compressed for semi-major axes between 200 and 400 pixels because the NUV and FUV colours change slowly for semi-major axes 200. The dashed lines show the UV colours for black-bodies and Castelli & Kurucz stellar atmosphere models of varying temperature along the line, with different extinction values for each line.

Simple Models for UV colours of the M31 bulge
The UV colour change of the bulge is shown in Figure 6. The semi-major axis, a, is encoded in a standard RGB colour scheme, with bluest for a = 0 gradually changing to green at a 50 , green changing to orange at a 100 and orange changing to red at the outermost analyzed ellipse, with a = 167 .
The larger errors and higher spread for the innermost ellipses (upper points in the diagram) are caused by the fewer pixels per annulus for annuli of small semimajor axis. The colour-colour diagram shows a significant change in the ultraviolet colour of M31's bulge. Points nearest the centre have a greater ratio of F148W flux to N219M flux than those further out, while having a similar ratio of N219M to N279N. Thus the FUV-NUV colour (vertical axis) of the bulge becomes "bluer" towards the core. Most of the colour change occurs between the centre and half-way out the analyzed region, from semi-major axis of 0 to 83 .
Colours for blackbody models and Castelli-Kurucz stellar atmosphere models (Castelli & Kurucz 2003, from the STScI website https://www.stsci.edu/ ) are plotted for comparison. For the stellar atmospheres, surface gravities were taken as the recommended values for each temperature, and solar metallicity was used. Model colours are calculated as the ratio of filter fluxes determined from a given model spectrum multiplied by the applicable UVIT filter response curve. Interstellar extinction is calculated using equations from Fitzpatrick & Massa (2007) which include the wavelength dependence for optical and UV. The extinction is parametrized by A V , given in magnitudes for the V band, and includes internal extinction within the bulge of M31 plus foreground extinction in the Milky Way along the line-ofsight toward M31.
The data are consistent with blackbodies between 9,500 and 11,000 K, with extinction values, A V , between 0.4 and 0.7. The stellar atmosphere models suggest a lower temperature range, 8500 to 9000 K, and a higher extinction, between 0.7 and 1.1. The stellar models are more realistic if the emission is from stellar atmospheres, particularly if the UV emission primarily comes from a population of unresolved hot stars. This change in model temperature can be interpreted as a rise in the average temperature stars towards the centre. Leahy et al. (2018) found evidence for a young hot stellar population in M31's bulge, and this combined with other evolved stars are possible causes of the colour changes.

Stellar Population Models for UV magnitudes of the M31 bulge
To construct more realistic models for the FUV and NUV flux from the M31 bulge, which consists of large numbers of unresolved stars, it is necessary to consider stellar population models. The basic models are simple stellar populations (SSPs) which consist of a set of stars of a single age and metallicity but with a range of initial masses. More complex models can then be constructed by combinations of these SSPs.
Here, we use the SSPs using the Padova stellar models, calculated using the CMD 3.4 online tool at http://stev.oapd.inaf.it/cgi-bin/cmd. In CMD 3.4, we used the recommended options: PARSEC evolutionary tracks (Bressan et al. 2012), version 1.2S, for premain sequence to first thermal pulsation or carbon ignition, and COLIBRI (Pastorelli et al. 2020) for thermal pulsation-asymptotic giant branch evolution, ending at total loss of the envelope. Circumstellar dust properties were used for M stars and for C stars from Groenewegen (2006), long-period variability along the red giant branch and asymptotic giant branch phases was taken from Trabucchi et al. (2019), and the initial mass function of Kroupa (2001) was used. For our application, synthetic photometry was calculated in AB magnitudes for the UVIT filter bands.
We add extinction, with parameter A V , using the extinction curve of Fitzpatrick & Massa (2007). The resulting SSP UVIT colours are compared to the colours  Table 2) of M31. T is the age in years and Z is the metallicity of the population. None of the SSPs can fit the data. Right panel: data for the Whole Bulge region (black symbols and line) and for the innermost annulus (0 to 58 in semi-major axis, blue symbols and line); model of 3 combined SSPs for the Whole Bulge region (red symbols and dashed line); model of 2 combined SSPs for the innermost annulus (magenta symbols and dashed line). Figure 9. Best-fit extinction (AV ) vs. outer semi-major axis (aout) of each annulus from the 2 SSP fits (red points) and from the 3 SSP fits (black points). Table 2 gives the inner and outer semi-major axes of each annulus, ain and aout. The two points for largest aout are for the outer-most annulus (lower point) and the Whole Bulge region (upper point). Linear functions are fit to AV vs. aout for all annuli for the 2 SSP model and fit to all but the innermost annulus for the 3 SSP model. of the M31 bulge in Figure 7. It is seen that the inner bulge colours (blue points) are similar to SSPs with A V ∼0.4-0.6 and a range of ages or metallicities. The middle bulge colours (orange points) are similar to SSPs with A V ∼0.2 and low age (log(T)∼ 6.5 to 7 and low metallicities (∼ −1.5 to -2).
The above comparison of SSPs with colours is incomplete, so we carry out fitting of the M31 bulge FUV-NUV spectral energy distribution (SED) to multiple stellar populations. This has the advantages of including data from all five UVIT filters, instead of three, and preserves the normalization which is required to obtain the size of the stellar populations. We extracted the FUV and NUV magnitudes for the whole bulge, here taken as an outermost ellipse with semi-major axis of 450 pixels (188 ). This large ellipse was divided into 10 annuli of equal-area to examine spatial variations. The semi-minor axis, b, to semi-major axis, a, ratios (b/a) were taken from the above analysis. For the whole and 10 annuli regions, systematic errors in conversion from counts/s to flux and magnitude are comparable to the 1/ √ N statistical errors, so were included. Table 2 shows the parameters of the elliptical annuli. For each area the UVIT magnitudes were fit by minimizing χ 2 with either 1, 2 or 3 SSPs. The fitting was carried out for each SSP in the grid of log(T) and log(Z/Z ). For each SSP the free parameters were mass of the SSP and A V . For the case of 2 or 3 SSP fits, the extinction was taken to be the same for both (or three) SSPs. The best-fit was chosen to be the SSP with the lowest χ 2 of all SSPs in the grid. For the case of 2 SSPs and 3 SSPs the grids were 4 and 6 dimensional, resp.
None of the single SSP models are able to fit the data. The left panel of Figure 8 shows a few of the single SSPs Figure 10. Left panel: Best-fit SSP masses vs. outer semi-major axis (aout) of each annulus from the 2 SSP fits, with red points for the older metal-poor population and blue points for the younger metal-rich population. Table 2 gives the inner and outer semi-major axes, ain and aout, of the 10 annuli and of the Whole bulge region. Right panel: Best-fit SSP masses vs. outer semi-major axis (aout) of each annulus from the 3 SSP fits, with red points for the oldest metal-poor population, cyan points for the younger metal-rich population, and blue points for the intermediate age and metallicity population. For each set, there are two points at largest aout: the lower mass is for the outermost annulus and the higher mass is for the Whole Bulge region. Exponential functions are fit to AV vs. aout for all annuli for the 2 SSP and 3 SSP models.
compared to the data for the Whole region (black line labelled M31 Bulge). The data errors are comparable to the symbol sizes. The SSP normalization was chosen so that the calculated magnitudes of the SSP did not exceed the data. It is seen that young SSPs (log(T) 9) are too faint for the longer wavelength data (by ∼2 magnitudes) and that old SSPs are too faint for the shorter wavelength data.
Including a second SSP improves the fit significantly. The right panel of Figure 8 shows the best-fit 2 SSP model (magenta line) for the innermost annulus/ellipse (data shown by blue symbols and solid line). Adding a third SSP further improves the fit. The Whole region data (top, black symbols and solid line) and the best-fit 3 SSP model (red symbols and dashed line) are shown. The best-fit masses and extinction (A V ) for the 2 and 3 SSP models are listed in Table 2.
For the 1 SSP model there are 4 parameters (2 grid parameters and 2 continuous parameters) and 5 data points. Best-fit χ 2 values are between 1000 to 2000 for the 10 annuli and whole regions and rule out single SSP models/ For the 2 SSP models there are 7 parameters (4 grid, 3 continuous), more than the number of data points, thus a statistically good fit should have χ 2 ∼ 1. The χ 2 values are between 4 and 9, indicating a large improvement over the 1 SSP model. The 3 SSP models have 10 parameters (6 grid, 4 continuous) and χ 2 ∼ 1 as expected for good fits. The 3 SSP models show significantly better fits than the 2 SSP models, but they have the disadvantage of having more parameters than the number of data points in the SED.

DISCUSSION
M31's bulge, analyzed here out to ∼ 170 , is elliptical. The ellipticity is different in the different NUV and FUV filters (Figure 3) with a smaller ellipticity found in the longer wavelength filters (e.g. at semi-major axis of 25 , a/b 1.17 in NUVN2 vs. 1.27 in F148W). The ellipticity is lowest near the centre: a/b = 1.14: this is apparent on the contour plots ( Figure 2). It increases outward, peaking at ∼ 17 major axis in the NUV, and at ∼ 25 in the FUV. Further from the centre the ellipticity decreases outward, to a minimum near ∼ 55 , then increases outward to the maximum semi-major axis analyzed here (∼ 95 ).
Outside of 17 , the ellipticity increases with decreasing wavelength (left panel of Figure 3). Inside of 17 the ellipticity decreases with decreasing wavelength, i.e. it is approximately reversed, likely indicating a different physical component, or different mix of stellar populations, in the central part of the bulge.
The bulge's inclination changes smoothly with increasing semi-major axis, increasing from 23 • reaching a peak of 48 • in all filters at ∼ 25 . It then slowly decreases outwards to 33 • , implying a clockwise rotation from centre to semi-major axis of 17 , followed by a gradual counterclockwise rotation from there outward to the maximum semi-major axis of 95 . This Note-a: The other best-fit parameters were log(T1)=10, log(Z1/Z )=-2, log(T2)=9, log(Z2/Z )=-1.5 for Ann1 to Ann10; for Whole the best-fit parameters were the same except log(Z2/Z )=-2.
Note-b: ain and aout are the inner and outer semi-major axes for each annulus.
small amount of rotation (25 • ) is in the contour plots ( Figure 2) but it is not easy to discern. The bulge is significantly more circular very near the centre. The brightness profile of the bulge vs. semi-major axis was extracted using 1 pixel wide elliptical annuli. There is a sharp increase of brightness towards the centre. The profile was fit with a Sersic function, showing an excess at the centre which is more extended ( 2 ) than a point source (right panel of Figure 4).
Thus we fit the bulge with the sum of a Sersic and a Gaussian and the sum of two Sersic functions. The double Sersic functions yield the best fit, although the decrease in χ 2 is highly significant only for the F148W data. We conclude that the profile is more complex than can be fit with simple analytic functions, such as we have used. For no galaxy other than M31 do we have such good photometry in NUV and FUV. Thus simple analytic fits to a galaxy bulge profile have not tested to this level previously.
The next step of analysis of the bulge considered spectral changes vs. semi-major axis distance. We analysed the colour-colour diagram using flux ratio of F148W to N219M (equivalent to N219M magnitude minus F148W magnitude) vs. flux ratio of N219M to N279N (equivalent to N279N magnitude minus N219M magnitude). In this diagram, hotter (brighter at shorter wavelengths) objects are higher in the diagram, and objects with more extinction are further left. The latter is a result of the peak of the UV extinction curve at 220 nm, matching the N219M filter.
From the colour-colour diagram (Figure 6), the centre of the bulge is seen to be hotter (blackbody temperature 11500 K, or Castelli-Kurucz stellar atmosphere temperature 8800 K) than the outer parts of the bulge (blackbody temperature 9750 K, or Castelli-Kurucz stellar atmosphere temperature 8450 K). The extinction of the bulge is between A V = 0.4 and 0.7 for blackbodies, or between A V = 0.7 and 1.1 for stellar atmospheres. The data errors are such that the outer bulge has distinctly lower extinction than the middle bulge, but the central bulge has extinction consistent with either middle or outer bulge.
Because the bulge NUV and FUV emission is from a large number (of order ∼ 10 9 ) of stars, we compare the colour-colour diagram to the colours calculated for single stellar populations (SSPs). SSPs with different ages and metallicities are shown. Figure 7 shows that comparison. The apparent trend is that the outer bulge (red points) is more consistent with a young (∼ 10 6.5−7 yr) metal poor (log(Z/Z ) ∼ −2) SSP, the middle bulge (green points) is more consistent with slightly less young (∼ 10 7−7.5 yr), less metal poor (log(Z/Z ) ∼ −1 SSP, and the inner bulge (blue points) is more consistent with older (∼ 10 8.5 yr), metal poor (log(Z/Z ) ∼ −1.5) SSP.
We carried out fits to the bulge spectral energy distribution (SED) using all five UVIT filters. No single SSP is a good fit (Figure 8). This shows that the colourcolour diagram does not contain enough information to properly constrain stellar populations. A mixture of two SSPs or three SSPs was able to fit the SED of the M31 bulge. The three SSP fit is significantly better than the two SSP fit.
The SEDs were extracted for 10 elliptical annuli in the bulge each with 10% of the whole bulge area, with outer semi-major axis of each annulus ranging from ∼ 60 to ∼ 190 . For each annulus two SSP and three SSP fits were carried out. The results are shown in Figure 9 for extinction, and Figure 10 for masses of the SSPs. There is a mixture of at least two and probably three different SSPs, based on χ 2 of the fits in Table 2. This holds for all 10 annuli, i.e., throughout the M31 bulge.
The extinction is variable, decreasing from centre (< 60 ) to outer part of the bulge (180 to 190 ) in a nearly linear manner. The slope from the two SSP fit is steeper than for the three SSP fit, but the extinction falls in a narrow range (A V = 0.4 to 0.5) for both models. The three SSP fits show evidence for a decrease in extinction for the centre.
The two SSP fit has an old (10 10 yr) metal poor (log(Z/Z ) = −2) population plus an intermediate age (10 9 yr) less metal poor (log(Z/Z ) = −1.5) population for the whole bulge and for all of the individual annuli. The masses of the old and intermediate age populations vary with distance from center but the ratio of intermediate SSSP mass to old SSP mass is typically 0.1. Both populations mass vs. distance from centre are well fit by exponential functions, with scale lengths of 65 for the intermediate age SSP and 110 for the old SSP. The higher concentration of the younger hotter population toward the centre can explain the trend seen in the colour-colour diagram.
The three SSP fit has an old (10 10 yr) metal poor (log(Z/Z ) = −2) population plus an intermediate age (10 9.5 yr) metal poor (log(Z/Z ) = −2) population and a young (10 8.5 yr) higher metallicity (log(Z/Z ) = −0.5) population. The masses of the three populations vary with distance from center. The ratio of intermediate SSP mass to old SSP mass is typically 0.15. The ratio of young SSP mass to old SSP mass is typically 0.02. The SSP masses vs. distance from centre are well fit by exponential functions, with scale lengths of 75 for young SSP, ∼ 450 for the intermediate SSP and 110 for the old SSP.
The dominant old metal poor population has the same scale height and normalization from the two SSP and three SSP fits. The three SSP fit (with significantly improved χ 2 ) effectively replaces the 10 9 yr component of the two SSP fit with two separate SSPs of ages 10 9.5 yr and 10 8.5 yr. The scale-heights of the latter two populations are very different: the 10 9.5 yr SSP has scaleheight 450 vs. 75 for the youngest SSP. The different scale-heights of the 3 SSPs explains the observed colour change of the bulge: it is hotter near centre where the youngest SSP has the largest fractional contribution. Hammer et al. (2018) carried out hydrodynamical simulations of mergers with M31 and compared the simulations to several sets of observations of M31, including parameters of the Giant Stream, the M31 disc agedispersion relation, the long-lived 10 kpc ring, the recent and long-lived star formation event, and the slope of the halo profile. They found that a 4:1 major merger where the nuclei merged 2-3 Gyr ago, with first passage 7-10 Gyr ago, could explain many of the observed features of M31. McConnachie et al. (2018) analyzed data from the Pan-Andromeda Archaeological Survey of M31 to summarize the substructures in the stellar halo of M31. Their analysis shows the most distinctive substructures were produced by at least 5 different accretion events in the past 3-4 Gyr.
The minimum of 5 accretion events (McConnachie et al. 2018) can be consistent with merger with a single major galaxy (Hammer et al. 2018) because there were a few passages between first passage and final merger of the nuclei. The passage would have created major streams of stars and gas which would cause accretion events at several different times and with different orbital parameters.
Our SSP fits to the bulge measure star formation in the bulge, and are consistent with the above results. This star formation is caused by accretion of gas or disturbance of ambient gas to induce star formation 4 , thus related to merger and accretion events. The 10 9.5 yr population comprises 17%, and the 10 8.5 yr population comprises 0.3% of the mass of the old population (Table 2). If there are more than 3 SSPs in the bulge, the current analysis is unable to detect them given the data errors.
There have been previous studies of the stellar populations in the bulge of M31 (Dong et al. 2018, Saglia et al. 2018 and references therein). Our results are consistent with theirs, in that the bulk of the mass of stars are old (age ∼10 Gyr), but not consistent in that they find that the old bulge population is metal rich, with [Fe/H]∼0.3 dex. The main difference between the previous studies and this study is that previous studies have used optical observations (optical HST photometry, or Lick indices), while we use FUV and NUV band photometry. One possibility is that the current set of SSP models are not accurate enough in the FUV to determine metallicity reliably. We plan to update the current study when SSP models which are better verified in the FUV become available.
The spatial structure of the M31 bulge has been studied by Courteau et al. (2011), with primary analysis of the Spitzer IRAC 3.6 micron observations of M31. That study fitted the bulge ignoring the inner few arcseconds containing the nucleus of M31. The fits used a Sersic profile for radial profiles along the major and minor axes or azimuthally averaged radially profiles. They used two fitting methods: non-linear least squares and Monte Carlo Markov Chain. The results (their Table 2) were a Sersic index, n, between 1.7 and 2.4 and an effective radius, R e , between 0.5 and 1.1 kpc. Both n and R e had typical errors of 10% . In comparison, they fit the 2-D image using the GALFIT software (Peng et al. 2010) and obtained similar results: n between 1.9 and 2.0 and R e between 0.9 and 1.0 kpc. Our longest wavelength is 279 nm, far shorter than the 3.6 micron data analyzed by Courteau et al. (2011). For 279 nm we find n values of 2.33, 2.30 and 2.27 for single Sersic, Sersic plus Gaussian or double Sersic fits, respectively (second, third and fourth sections of Table 1). These values are only slightly higher than those from Courteau et al. (2011). For the shorter NUV 219 nm and the three FUV bands, we find significantly higher Sersic indices (3.2 to 5.6). The fact that the Sersic indices are quite different for NUV and FUV should not be too surprising, given that the shorter wavelengths are highly sensitive to the rare young and hot stars. These we have shown are more concentrated in the centre of the bulge (most clearly shown here in Fig. 10).

CONCLUSIONS
We have carried out observation of the bulge of the Andromeda Galaxy in NUV and FUV wavelengths using the UVIT instrument on the AstroSat Observatory. M31's bulge is elliptical in shape with small but distinctive changes in ellipticity and position angle with distance from the nucleus of M31.
The bulge is centrally concentrated in UV, and approximately follows a Sersic profile. The nucleus (in the innermost 4 ) is complex, as known from optical studies (Lockhart 2017). For the brightness profiles in all filters but F148W, the bulge well fit by Sersic plus Gaussian or double Sersic functions. For F148W, with its high signal-to-noise, none of the attempted functions gave statistically good fits, but the double Sersic had the best-fit, and gave no systematic residuals vs. radius, which may indicate that the bulge has small scale structure or clumps which are not fit by any smooth function.
The UV colour-colour diagram shows that the bulge becomes significantly bluer towards the centre. We carried out an analysis of the FUV -NUV SEDs using SSPs for the whole bulge and the bulge subdivided into 10 elliptical annuli. Our best fit model has three SSPs at each radius and yields the same best-fit ages and metallicities vs. radius for these SSPs. There is a dominant old SSP (10 10 yr, log(Z/Z ) = −2), an intermediate-mass, intermediate-age SSP (10 9.5 yr, log(Z/Z ) = −2) and a young low-mass SSP (10 8.5 yr, log(Z/Z ) = −0.5). The different SSPs have exponential mass profiles vs. radius with different scale-heights for each SSP. The ages and masses of the SSPs are consistent with an active merger history for M31, as known from studies of M31 using optical data.
The FUV and NUV data from UVIT/AstroSat are useful to measure the young stellar populations which emit at these wavelengths. Future studies of regions beyond the bulge promise to turn up more important properties of the ages and spatial distributions of stellar populations in M31.