© ESO 2020

1 Introduction

The evolution of spiral galaxies is driven by star formation and the matter cycle between the stars and the interstellar medium (ISM). Massive OB stars inject energy into the ISM through their radiation, stellar winds, and, finally, by supernova (SN) explosions. These processes are often correlated in space and time, producing superbubbles with sizes of typically 100–1000 pc. However, as we are located inside the Galactic disc, we are not able to observe the ISM in the Milky Way with all these processes and resulting structures in their entirety. Instead, for a detailed study of the ISM, we have to look beyond our Galaxy.

The Andromeda galaxy (M31) is the largest galaxy in the Local Group and the nearest spiral galaxy to the Milky Way, located at a distance of 783 kpc (Conn et al. 2016). With a size and a mass comparable to that of our Galaxy, this archetypal spiral galaxy provides a unique opportunity to study and understand the nature and the evolution of a galaxy similar to our own. Based on Spitzer data, Gordon et al. (2006) showed thatM31 has spiral-arm structures merged with the prominent star-forming ring at a radius of ~ 10 kpc. The Herschel images show a radial gradient in the gas-to-dust ratio and indicate that there are two distinct regions in M31 with different dust properties inside and outside a radius R ≈ 3.1 kpc (Smith et al. 2012; Fritz et al. 2012). Galvin et al. (2012) and Galvin & Filipovic (2014) presented high resolution radio continuum images of M31 from the Very Large Array, identifying 916 unique discrete radio sources across the field of M31 and faint extended emission at the location of the 10 kpc ring.

The star-formation history (SFH) in M31 has been studied in detail in observations with both the Hubble Space Telescope (HST, Lewis et al. 2015) and large ground-based telescopes (e.g. the Local Group Galaxy Survey, LGGS, Massey et al. 2006). Williams (2003) measured a mean star-formation rate of about 1 M_⊙ yr⁻¹ in the full disc of M31 over the last 60 Myr and produced maps of star-formation rate in different age ranges. Furthermore, the northern disc was observed in the Panchromatic Hubble Andromeda Treasury (PHAT) survey (Dalcanton et al. 2012) resolving over 100 million stars in M31 in the near-infrared (NIR) to the ultraviolet (UV). Using these data, Lewis et al. (2015) showed that the SFH varies significantly on small spatial scales in the northern disc. They found that the 10 kpc ring is at least 400 Myr old with ongoing star formation and estimate about 60% of all star formation in the last 500 Myr occuring within the ring. By extrapolating the results from the PHAT survey area, the total star formation rate in the disc was estimated to be ~ 0.7 M_⊙ yr⁻¹ over the last100 Myr, consistent with previous broadband estimates.

In X-rays, a survey of the entire M31 was performed with XMM-Newton (Stiele et al. 2011). These XMM-Newton observations have revealed an extended diffuse emission in M31, in particular on the northern side of the galaxy, which is well correlated with the star-forming ring. We carried out additional deep observations in two fields of the northern disc to perform a detailed study of the X-ray source population and the diffuse emission. Based on the new XMM-Newton data, we created a catalogue of X-ray sources down to luminosities of ~ 7 × 10³⁴ erg s⁻¹ (0.5–2.0 keV) and improved the classification of the previously known X-ray sources. The results of the point sources study (catalogue, spectral and timing studies, multi-wavelength cross-correlation, etc.), based on the new XMM-Newton observations, have been published by Sasaki et al. (2018, hereafter Paper I), in combination with the results of surveys performed with the Chandra X-ray Observatory and the PHAT data taken with the HST (Williams et al. 2018). In this paper, we present an analysis of these new deep XMM-Newton pointings (~ 200 ks on each field) focusing on the diffuse emission in order to study the structure and the physical properties of the hot ISM plasma. The interstellar space in galaxies like the M31 and the Milky Way are filled with cool clouds (T ≲ 10² K) of neutral hydrogen embedded in warm (T ≈ 10⁴ K) intercloud medium of partially ionised hydrogen (Cox 2005, and references therein). At the distance of M31, 1′′ corresponds to 3.8 pc, so the spatial resolution of XMM-Newton allows us to resolve the hot gas in bubbles and superbubbles in M31 on 100–1000 pc scales and study different morphologies, which might be indicativ of inhomogeneities in the ISM. The spectral resolution of the EPIC instruments allows us to investigate variations in the physical properties of the ISM such as luminosity, temperature, and metallicity, as well as absorption by cold ISM in the foreground. The distribution and the filling factor of the hot phase of the ISM can be studied in detail by comparing the X-ray emission to H I (e.g. Braun et al. 2009), Hα (LGGS) data, and cold gas and dust (Fritz et al. 2012).

We present the XMM-Newton data used for this work and the analysis methods applied to the study of the diffuse X-ray emission in Sect. 2. Also, we discuss the contribution of the unresolved point sources to the observed diffuse X-ray emission. In Sect. 3, we give an overview of the additional multi-wavelength data. In Sect. 4, morphological studies and spectral analyses are presented, followed by a study of the stellar population in each of the regions and the energy budget in Sect. 5. A summary and conclusions are presented in Sect. 6.

2 Observations and data reduction

We performed deep observations of two fields in the northern disc of M31 with XMM-Newton (Jansen et al. 2001) as part of our Large Project survey (PI: M. Sasaki), details of which can be found in Paper I (see also Table A.1). In the following sections we describe the reduction of these data and supplementary X-ray data for our imaging and spectroscopic analyses.

2.1 Imaging

In addition to our new deep fields, XMM-Newton has observed M31, including the northern disc, on many occasions over its lifetime. The majority of these observations focus on the core region, but a survey of M31 has also been performed covering the ellipse enclosed by the D₂₅ isophote down to a limiting luminosity of 10³⁵ erg s⁻¹ in the 0.2–4.5 keV band (Stiele et al. 2011). In addition, various interesting objects located around M31 have also been covered in separate, dedicated observations. In order to properly reveal and characterise the hot ISM plasma in M31 we made use of all available XMM-Newton data taken up to and including our deep exposures of the northern disc. The area of XMM-Newton coverage used for our M31 mosaic corresponds to that of Stiele et al. (2011). This coverage is plotted in relation to the stellar population in their Fig. 1, with regions N1 and N2 roughly corresponding to our new deep fields in the northern disc.

To produce mosaics of M31 and its northern disc we used the mosaicing capabilities of the XMM-Newton Extended Source Analysis Software (XMM-ESAS), packaged in SAS¹ 15.0.0. XMM-ESAS is based on the software used for the background modelling described by Snowden et al. (2004), and comprises a set of tasks to produce images from observational data and to create model quiescent particle background (QPB) images that can be subtracted from the observational science products (see Kuntz & Snowden 2008; Snowden et al. 2008).

We considered only those observations where the primary instrument was the European Photon imaging Cameras (EPICs), which use a pn-type CCD (Strüder et al. 2001) or a MOS-type CCD (Turner et al. 2001) imaging spectrometers (one EPIC with pn CCD and two with MOS CCD). We used only EPIC observations taken in Full Frame mode and processed each of the observational datasets according to the ESAS Cookbook². Standard filtering and calibration were applied using the epchain, emchain, pn-filter, and mos-filter tools. Finally, we discarded EPIC-pn and EPIC-MOS observations with soft proton (SP) flare-filtered exposure times less than 15 ks. Table A.1 summarises all the data included in our final screened sample.

The XMM-ESAS task cheese-bands was used to produce source exclusion masks for each observation. This task performs the source detection on all available EPIC instruments in multiple energy bands for a given observation and determines source exclusion regions accounting for the shape of the point spread function at the source location on the detector. We chose the 0.3–2.0 keV and 2.0–7.0 keV bands so that the source detection would be sensitive to both soft and hard sources and defined the exclusion limit as 1.1 times the W90 radius (the 90% enclosed energy radius). Source lists from individual observationswere merged using the merge_source_list task, and a final source mask was created using the make_mask_merge task. Bright sources at or near the edges of individual fields-of-view (FOVs) presented a problem, however. Uncertainty in source centroids resulted in the merge_source_list task not recognising that what appeared to be two separate sources were in fact the same source. In these cases we manually added the exclusion regions to the masks.

The pn-spectra and mos-spectra tasks were then used to produce images in the 0.3–0.7 keV, 0.7–1.25 keV, and 2.0–4.5 keV, with point sources masked. We excluded the 1.25–2.0 keV range to remove the contribution of the strong Al Kα and Si Kα instrumental fluorescence lines from the images. The pn-back and mos-back tasks were then used to produce the corresponding QPB images. Using XMM-ESAS, it is also possible to model and subtract residual soft-proton (SP) contamination from the images. This is important when producing large mosaics as some observations are affected by residual SP while others are not. To this end we extracted full-field spectra and QPB backgrounds for each observation and fitted the spectra in the 3–7 keV energy range where the SP contamination is expected to dominate. In observations where residual SP contamination was identified, the determined spectral parameters were used to generate model SP contamination images with the XMM-ESAS task proton. Subsequently, we used the merge_comp_xmm task to create count, exposure, QPB, and residual SP mosaic images. Finally, the adapt_merge task was used to create particle background subtracted and exposure corrected mosaics in each energy band, binned into ~ 10′′ × 10′′ bins and adaptively smoothed to yield the final images. To aid in the visualisation of the diffuse emission in M31, “holes” due to masked point sources were filled using the biharmonic function based inpainting algorithm³ contained in the scikit-image python package⁴. These mosaics are presented in Fig. 1.

Fig. 1 XMM-Newton mosaic images of M31 with 0.4–0.7 keV shown in top-left and 0.7–1.25 keV in top-right. The white dashed circles mark the fields of our northern disc LP survey, which are revealed in more detail in the zoomed images in bottom-left and bottom-right. The XMM-Newton coverage in relation to the stellar population is shown in Fig. 1 of Stiele et al. (2011). The colours are in linear scale as shown by the colour bars in each image. The units of the images are cts s−1 deg−2.

Fig. 2 Left: XMM-Newton 0.4–1.25 keV image of the northern disc of M31 with contours. The image has been adaptively smoothed and masked point sources filled as described in the text. Yellow, cyan, blue, and red contours correspond to 5, 10, 20, and 50σ above the average background level determined from regions exterior to the northern disc. Regions of interest discussed in the text are marked. Right: XMM-Newton 0.4–1.25 keV image of the northern disc with the spectral extraction regions defined from the contours in left shown in cyan. Background regions are shown in yellow. The deep EPIC northern disc fields are shown by the green circles.

2.2 X-ray morphology

To identify regions of M31 with significant X-ray emission, we determined the average background surface brightness in our mosaic and corresponding standard deviation (σ) in the 0.4–1.25 keV band from several background regions exterior to the northern disc. We then defined contours at 5σ, 10σ, 20σ, and 50σ above the background level. The resulting image and contours are shown in Fig. 2-left. The X-ray contours follow the spiral structure of M31 (see discussion on multi-wavelength morphology in Sect. 5.1 and Fig. 8).

Three bright regions of X-ray emission are evident along the bright star-forming ring which are most likely related to high mass stellar populations (see Sects. 5.1 and 5.2.5 for further details). We labelled these as Regions 1–3 in Fig. 2-left. A fourth bright emission region is seen closer to the core of M31 though still appearing to be associated with a massive stellar population. We labelled this as Region 4 (Fig. 2-left). Larger, fainter extended emission is observed just inside the star forming ring, filling the gap between the ring and inner dust lane (see Fig. 8). This emission is broken towards the north. Therefore, we separated this ring-adjacent emission into two called Regions 5 and 6, as shown in Fig. 2-left. The faint extended emission inside the dust lanes was labelled Region 7. We used the 5σ contours to define spectral extraction regions for these extended emission regions, which are shown in Fig. 2-right.

The brightest emission in our mosaic was found in the nuclear region. While some of this emission is likely diffuse in origin (Li & Wang 2007), much of the emission is due to the unresolved population of X-ray emitting objects such as cataclysmic variables (CVs), active binaries (ABs), X-ray binaries (XRBs), etc. In addition, the multitude of very bright and variable X-ray sources in and around the nucleus, and the relatively large XMM-Newton point-spread-function, makes the spectral extraction and analysis of this region very challenging. For this reason, we do not analyse the emission from the M31 core.

2.3 Spectra

The mosaics constructed from the available data showed diffuse structure in the northern disc of M31 only. This is due to our new deep observations of this region and highlights the necessity of similarly deep observations in the southern disc to reveal its hot ISM plasma. Therefore, our spectral analysis is focused solely on the northern disc viewed through our observations. We note that the spectral analysis presented in this work is somewhat complex with multiple thermal and non-thermalcomponents contributing to a single spectrum with these having instrumental, particle, and astrophysical origins. Therefore, particularly given the faint nature of the emission being studied, it is essential to constrain the various contributions to our spectra as much as possible. In this section we outline our approach to extracting and analysing the spectra, and constraining the instrumental, astrophysical background, and unresolved M31 source contributions.

For the spectral extraction we closely followed the procedures discussed in detail in Maggi et al. (2016), only differing in our approach to the background fitting. In this work, we first constrained the background emission, then fixed these emission components in fits of the extended emission region spectra, rather than fit spectra of extended and background regions simultaneously. We used standard tools in SAS to produce our spectra and associated files. Prior to extracting spectra, we generated vignetting-weighted event lists for the EPIC instruments to correct for the effective area variation across the FOV using the SAS task evigweight. Spectra were extracted from our extended emission regions (see Sect. 2.2), with background spectra extracted from regions within our deep fields (see Sect. 2.3.1). Out-of-time spectra were also created, scaled, and subtracted from the extended and background spectra following standard SAS procedures⁵. The point source catalogue from Paper I was used as input to the SAS task region to determine source exclusion regions. Exclusion regions were defined using the 90% encircled energy fraction of the sources. Extended and background spectra were extracted using the SAS task evselect, with redistribution matrices and ancillary response files produced by the SAS tasks rmfgen and arfgen, respectively. The spectra were grouped to minimum of 30 counts per bin to allow the use of the χ² -statistic. All fits were performed using XSPEC (Arnaud 1996) version 12.9.1 with ATOMDB version 3.0.3⁶. Abundance tables were set to those of Wilms et al. (2000) and photoelectric absorption cross-sections to those of Balucinska-Church & McCammon (1992).

2.3.1 X-ray background

We were somewhat limited in our choice of background regions in the northern disc fields given that our extended emission regions cover most of the EPIC instrument FOVs. We were further constrained as we required regions that fell on the detectors in both deep pointings at each northern disc field (see Fig. 2-right green circles). This is important for identifying any time variable background components, such as residual SP contamination and Solar Wind Charge Exchange (SWCX) emission. Three background regions were chosen which we hereafter refer to as SE, NE, and NW, and are shown in Fig. 2-right. Spectral fits were performed on these regions to characterise the X-ray background components.

Astrophysical components

The astrophysical X-ray photon background (AXB) typically comprises four or fewer components (Snowden et al. 2008; Kuntz & Snowden 2010), namely the unabsorbed thermal emission from the Local Hot Bubble (LHB, kT ~ 0.1 keV, fixed in the model), an absorbed cool Galactic halo (kT_cool ~ 0.1 keV), a higher temperature (kT_hot ~ 0.2−0.7 keV), absorbed thermal component representing emission from the hotter Galactic halo or intergalactic medium (IGM), and an absorbed power law representing unresolved background active galactic nuclei (AGN). To account for the unresolved AGN in our fits, the photon index (Γ) and normalisation of the power law were fixed to Γ = 1.46 and the equivalent of 10.5 photons cm⁻² s⁻¹ sr⁻¹ at 1 keV (Chen et al. 1997), respectively. In the case of our observations, halo emission from M31 may also supplement the Galactic halo and IGM emission. Rather than explicitly model this, we simply let the two thermal components represent the combined emission from the Galaxy, M31, and the IGM. All thermal components were fitted with the apec (Smith et al. 2001) thermal plasma model in XSPEC. Photoelectric absorption models (the model in XSPEC) were included to account for absorption in both the Galaxy and M31. The Galactic value was fixed to 1 × 10²¹ cm⁻² based on theDickey & Lockman (1990) H I maps, determined using the HEASARC N_H Tool. We also note that X-ray photons resulting from SWCX can affect the fitting and interpretation of thermal plasma model applied to sources of interest. This time variable emission results from solar wind ions interacting with neutral atoms in the heliosphere or Earth’s magnetosheath. This emission is characterised by emission lines from the ion species in the solar wind, such as C, O, Mg, etc., which can masquerade as thermal emission lines. In trial fits to our spectra we searched forany variation between pointings at our northern disc fields but found no evidence of emission line variation. We concluded that SWCX does not significantly affect our spectra and ignored it during subsequent analyses.

Particle induced components

The particle-induced background of the EPIC consists of the QPB, instrumental fluorescence lines, electronic read-out noise, and residual SP contamination. To determine the contribution of the first three components we made use of vignetting corrected filter wheel closed (FWC) data. Due to the very strong and position dependent contribution of electronic noise below 0.4 keV and the strong fluorescence emission line contribution above 7 keV for the EPIC-pn detector, we limit our spectral analysis to the 0.4–7 keV range. We extracted spectra from FWC data from corresponding regions on the detector as our extended regions and the backgrounds. The EPIC-pn and EPIC-MOS FWC spectra were fitted with the empirical models developed by Sturm (2012) and Maggi et al. (2016), respectively. Since these spectral components are not subject to the instrumental response, we used a diagonal response in XSPEC. The FWC data were fitted simultaneously with the corresponding source and background spectra, with FWC component parameters linked between the two. The only variables were the widths and normalisations of the strongest fluorescence lines, as well as a scaling factor to account for any difference in the QPB level between the FWC and extended and background spectra. Residual SP contamination is a time and position dependent component of the particle background. We allowed for this by including a power law not convolved with the instrumental response in the fits to the source and background spectra Kuntz & Snowden (2008).

We fitted the background spectra simultaneously with their corresponding FWC spectra to account for the particle-induced counts. Multiplicative factors accounting for the area of the spectral extraction regions in arcmin² were included in the XSPEC model, normalising the astrophysical model components to acrmin². These factors were derived from the BACKSCAL keyword set in the metadata during processing, which takes into account missing sky area due to chip gaps, bad pixels, etc. The fit results are given in Table 1 with the fits shown in Fig. 3. We found that the spectra between the three background regions were consistent with each other. Even though we included a 0.1 keV thermal plasma model to account for any LHB contribution, this component tended to 0 in each region. This is due to our restricted energy range of 0.4–7 keV, making our fits somewhat insensitive to this component. The defined background regions are ideal for some of our extended emission regions, but not all. In particular, Regions 4, 6, and 7 are not located adjacent to any (see Fig. 2). This is because of our background selection criteria described above. To account for the background in these regions and given that the NE, NW, and SE backgrounds are consistent with each other, we generated a master background by fitting the spectra of all background regions simultaneously in XSPEC, the results of which are given in Table 1 (Master) with the fits shown in Fig. 3, bottom-right. These best fit parameters were used to fix the astrophysical background contribution in the extended Regions 4, 6, and 7.

Fig. 3 XMM-Newton spectral fits of the background regions. Top-left: simultaneous spectral fits to the background region NE, and corresponding FWC spectra extracted from EPIC-MOS2 (green and blue) and EPIC-pn (black and red). The dashed lines show the combined AXB contribution while the dotted lines show the combined particle-induced contribution to the extended emission spectra. For clarity, only spectra from ObsID 0763120301 are shown. Top-right: same as top-left but for background region NW. Bottom-left: same as top-left but for background region SE. In this case the spectra from ObsID 0763120101 are shown. Bottom-right: simultaneous fits to all background spectra with each background region source and FWC spectrum shown in a different colour. The combined AXB and FWC contributions are also shown for each. See text for details.

Fig. 4 X-ray flux distributions of M31 sources detected in the southern (top row) and northern (bottom row) deep fields in the northern disc. Left panels: X-ray flux functions, along with power law and cut-off power law fits. Right panels: combined flux distribution expected from sources in flux bins, normalised for the deep field areas. The unresolved flux per unit area is determined by integrating over this distribution and subtracting the flux of detected sources.

2.3.2 Unresolved sources

In addition to the background and extended emission contributions, the extended region spectra will be contaminated by emission from unresolved point sources. This will comprise a multitude of X-ray sources which are intrinsically weak such as CVs and ABs, or the faint end of the XRBs, supernova remnant (SNR), and supersoft source (SSS) populations. It is important to account for these components in the spectral analysis to isolate the truly diffuse emission contribution. In the following we describe how we account for these in the spectral fits.

For the unresolved, though intrinsically bright X-ray sources, such as the XRBs, SNRs, and SSSs, we estimated the unresolved flux by constructing X-ray flux functions using the source catalogues of Paper I for each of the deep fields in the northern disc (see Fig. 4). Only sources with confirmed source types were included. To estimate the unresolved flux we first fitted the X-ray flux functions with either a power law or cut-off power law model (see Fig. 4, left panel of both rows). The cut-off power law provides a slightly better fit to sources in the southern field of the northern disc, closer to the nuclear region of M31. We therefore adopted this model to describe the northern sources, though the power law fit is included in Fig. 4 for completeness. For the northern field, the cut-off power law fit was indiscernible from the normal power law so we only consider the latter. Using these fits we then determined the combined flux per unit area from sources in each flux bin, normalised for the deep field area (Fig. 4, right panels). We integrated over this distribution to determine the total flux per unit area expected from the entire source population and subtracted the resolved component to obtain the unresolved flux per unit area. We found the unresolved populations for the southern and northern deep fields to be 1.8 ×10⁻¹⁵ erg cm⁻² s⁻¹ arcmin⁻² and 1.9 ×10⁻¹⁵ erg cm⁻² s⁻¹ arcmin⁻², respectively.These values were then used to estimate the normalisations for the unresolved model component in spectral fits to the extended regions. However, we also require a spectral shape for this component.

To estimate the spectral shape, we extracted spectra, backgrounds and response files from each of the sources used to construct the X-ray flux functions in Fig. 4 and combined them using the SAS task epicspeccombine. Source regions werethe same as the source exclusion regions described in Sect. 2.3.1 and background regions weredefined using annuli with inner and outer radii of 1.1 and 2.0 times the source extraction radius, respectively.We omitted very bright sources to prevent them from dominating the combined spectrum, setting a source cutoff of 500 counts per source. The combined faint source spectrum is shown in Fig. 5. We fitted the spectrum with several spectral models but found that a simple absorbed power law was sufficient with N_H = 0.96 (0.74−1.2) × 10²¹ cm⁻² and Γ = 2.37 (2.22−2.53) giving a reduced χ² = 1.07 for 662 degrees of freedom. Based on the best fit spectrum and expected contribution of the unresolved sources, we determined model normalisations for regions in the southern and northern fields of K = 5.86 × 10⁻⁷ photons keV⁻¹ cm² s⁻¹ arcmin⁻² and K = 5.55 × 10⁻⁷ photons keV⁻¹ cm² s⁻¹ arcmin⁻² at 1 keV, respectively. These spectral parameters were frozen in the fits to the extended region spectra.

We cannot apply the same method to determine the contribution from the intrinsically weak source population since none were detected in Paper I. These sources are associated with the older stellar population in M31. Revnivtsev et al. (2007, 2008) demonstrated that the X-ray luminosity per unit stellar mass due to weak X-ray sources (CVs and ABs) associated with the old stellar population in the solar vicinity (Sazonov et al. 2006) is compatible with that observed in other galaxies. Therefore, we can determine the expected X-ray luminosity from the unresolved CVs and ABs in the > 1 Gyr population using the mass or, equivalently, luminosity of the stellar population, which can be estimated using NIR data, specifically the K_s -band. We determined the combined K_s-band luminosity of our extended regions using the 2MASS Large Galaxy Atlas (Jarrett et al. 2003)⁷. We first identified and excised foreground stars from the images. We initially obtained estimates for each extended region. However, it became clear during the spectral analysis that this contribution was most significant in Regions 5, 6, and 7, but remained small when compared to both the background and diffuse components. Nevertheless, we included a model to account for unresovled CVs and ABs in our models for these three regions. We determined K_s = 9.6, 7.1, and 4.7 for Regions 5, 6, and 7, respectively. Converting these into a luminosities yields $L_{K_s}=1.7\times {10}^7$ L_⊙ for Region 5, $L_{K_s}=1.9\times {10}^8$ L_⊙ for Region 6, and $L_{K_s}=1.6\times {10}^9$ L_⊙ for Region 7. We then used the relation $L_{X,0.5-2\text{keV}}/L_{K_s}=(5.9\pm 2.5)\times {10}^{27}$ erg s⁻¹ $L_{K_s,L_{\odot }}^{-1}$ (Revnivtsev et al. 2008) and adjusted for the X-ray energy range using their emission model, a mekal model with kT = 0.5 keV and solar abundance, and a power law with Γ = 1.9, with the ratio of component luminosities equal to 2.03, to determine an expected contribution to the X-ray emission of L_{X,0.3−10 keV} = (1.1 ± 0.4) × 10³⁵ erg s⁻¹ for Region 5, L_{X,0.3−10 keV} = (1.1 ± 0.5) × 10³⁶ erg s⁻¹ for Region 6, and L_{X,0.3−10 keV} = (9.8 ± 0.4) × 10³⁶ erg s⁻¹ for Region 7. We then converted these values to flux and normalised for area. The resulting values were used with the mekal+powerlaw model of Revnivtsev et al. (2008) to determine the model normalisations K_mekal = 7.61 × 10⁻¹⁰ cm⁻⁵ arcmin⁻², K_mekal = 3.40 × 10⁻⁸ cm⁻⁵ arcmin⁻², and K_mekal = 9.16 × 10⁻⁸ cm⁻⁵ arcmin⁻² for Regions 5, 6, and 7, respectively. The power law normalisations were tied to the mekal normalisation so that their luminosity ratios were 2.03, as per Revnivtsev et al. (2008).

Overall, the unresolved point sources were found to contribute significantly to the observed flux from the extended emission regions. From the spectral fits described in Sect. 4, we estimated that point sources account for ~ 7% of the flux from Region 1, ~ 14% from Region 2, ~ 34% from Region 3, ~ 35% from Region 4, ~ 42% from Region 5, ~ 59% from Region 6, and ~ 71% from Region 7. The intrinsically weak sources contributed <1%, ~ 2%, and ~ 7% of the flux to Regions 5, 6, and 7, respectively.

Fig. 5 Combined faint source spectrum, fitted with an absorbed power law model with best fit NH = 0.96 (0.74−1.2) × 1021 cm−2 and Γ = 2.37 (2.22−2.53) giving a reduced χ2 = 1.07 for 662 degrees of freedom.

2.3.3 M31 emission

We expect any emission from the hot ISM plasma in M31 to be dominated by thermal processes. Therefore, we included thermal plasma models in the spectral fits absorbed by foreground Galactic and M31 material on top of the background components discussed in the preceding subsections. To fit the emission we proceeded through the following steps: we first tried a single-temperature thermal plasma in collisional ionisation equilibrium (CIE), namely the apec model in XSPEC; if a good fit to the data was found we tried a non-equilibrium ionisation (NEI) model (nei in XSPEC, Borkowski et al. 2001) to provide constraints on the ionisation timescale (τ)⁸ of the plasma; if the apec fit was bad, we tried the nei model to test for a non-equilibrium plasma and two temperature models (combinations of apec and nei) in case the emission comprised two distinct thermal components, though the latter was never required. The extended emission being studied is very faint in nature and, even with the deep observationsof the northern disc, the spectral statistics did not allow for any more complex models with variable abundance.

2.3.4 Neighbouring region backgrounds

It is clear from Fig. 2 (right) that there is some overlap of the larger extended emission regions with the more localised emission association with the young stellar populations, particularly along the dust ring. Namely, the extended emission in Region 5 extends across and likely contributes to emission in Regions 2 and 3. The extended emission in Region 7 is similarly affected by Region 4. To account for this we included a spectral component for Region 5 in fits to Regions 2 and 3, and from Region 7 in fits to Region 4. The contributions of these components were fixed to the best fits in Regions 5 and 7 and normalised per unit area.

3 Ancillary data

To compare the morphology and properties of the hot ISM plasma in M31 to the stellar and colder ISM components of the galaxy, we made use of multi-wavelength archival data presented in the literature, which we describe below.

3.1 UV

The Galaxy Evolution Explorer (GALEX) satellite has surveyed M31, obtaining images in the near-UV (NUV) and far-UV (FUV), the analysis of which is presented in Thilker et al. (2005). The angular resolution of GALEX allows for the separation of individual young stellar clusters hosting massive stars from each other, their environments, and the diffuse UV emission in the galaxy. Therefore, these observations provide a tracer of the populations of massive stars driving the dynamics of the ISM. We obtained the GALEX exposures from the Mikulski Archive for Space Telescopes⁹. We mosaicked the NUV and FUV fields using the Montage tool¹⁰.

3.2 Optical

M31 was observed as part of the Local Group Galaxy Survey (LGGS), performed using the Kitt Peak National Observatory 4 m telescope (Massey et al. 2006). The study covered 2.2 deg² along the major axis of M31 and includes all massive star forming regions. Broadband images were obtained in U, B, V, R, and I, as well as Hα, [S II ], and [O III] narrow band, emission line images. The warm ISM (10⁴ K) of M31 is revealed in Hα. Therefore, we downloaded both the bias-corrected and flat-fielded Hα and R-band images from the LGGS archive¹¹ in order to produce continuum subtracted Hα images to reveal the full extent of the faint diffuse warm ISM of M31. The continuum subtraction was performed on the individual fields of the LGGS survey (10 in total) before constructing a mosaic of these fields. The R-band image was scaled to the Hα image in each field by first identifying stellar sources in each using the DAOStarFinder task (based on DAOFind, Stetson 1987) contained in the Astropy (Astropy Collaboration 2013) affiliated package Photutils (Bradley et al. 2016), determining the scaling relation based on their fluxes, then scaling and subtracting the continuum image from the emission line image. The resulting continuum subtracted fields were mosaicked using the Montage tool.

3.3 H I

Neutral H in the M31 ISM can be traced using the 21 cm emission line. Therefore, we make use of the high-resolution H I cube presented in Braun et al. (2009), obtained using the Westerbork Synthesis Radio Telescope and the Green Bank Telescope. These data reveal the atomic gas component in M31 with a ~ 15′′ spatial and 2.3 km s⁻¹ velocity resolution in a 7 × 3.5 deg² region.

3.3.1 Infrared and sub-mm

The cold ISM of M31 can be revealed by infrared (IR) emission. We used data from the Spitzer Space Telescope (Werner et al. 2004) and the Herschel Space Observatory (Pilbratt et al. 2010). The 24 μm band of the Multiband Imaging Photometer (MIPS, Rieke et al. 2004) provides us with a picture of the stochastically and thermally heated dust in the ISM, and is presented in Gordon et al. (2006). We obtained these data from the IPAC Infrared Science Archive¹². Colder dust (~ 10 K) requires observations longward of 200 μm, in the wavelength range of the Herschel Spectral and Photometric Imaging Receiver (SPIRE, Griffin et al. 2010). Observations of M31 with Herschel SPIRE are reported as part of the Herschel Exploitation of Local Galaxy Andromeda (HELGA) survey by Fritz et al. (2012). For our study, we make use of the SPIRE 250 μm data to trace the cold dust morphology of M31.

3.3.2 Panchromatic Hubble Andromeda Treasury

The Panchromatic Hubble Andromeda Treasury (PHAT; Dalcanton et al. 2012) is a Hubble Space Telescope (HST) survey which observed one-third of the star-forming disc of M31 in six filters, from the UV to the NIR. The survey provided luminosities and colours for over 100 million individual stars. We made use of the stellar photometry collected in this survey (F475W and F814W filters, Williams et al. 2018) for the study of the stellar populations in our regions of interest.

4 Results

The spectral analyses of each individual extended emission regions are presented in the following sub-sections. For clarity, we only show spectral plots of the EPIC spectrum with the most counts for each, along with additive components of interest (Figs. 6 and 7). Plots showing the simultaneous fits to each observational and FWC spectrum, along with all additive components, are given in Appendix B for completeness (Figs. B.1 and B.2).

4.1 Region 1

We found that the emission over the background components in Region 1 could be fitted with a single apec model with kT = 0.20 (0.18−0.24) keV with N_H = 0.62 (0.39−0.76) × 10²² cm⁻² giving a reduced χ² (${\chi }_{\nu }^2)=0.99$. We also tried the nei model which provided a similarly good fit and found good agreement with the apec model, with the additional free parameter τ > 4.06 × 10¹¹ s cm⁻³ suggesting a plasma close to or in CIE. The fit results are given in Table 2, with the spectra and nei fits shown in Fig. 6, top left.

4.2 Region 2

We fitted the extended emission with a single apec model and found kT = 0.60 (0.56−0.65) keV with N_H = 0.07 (0.02−0.10) × 10²² cm⁻² giving a ${\chi }_{\nu }^2=0.99$ (see Table 2). We also tried the nei model, however the τ parameter was not well constrained rendering the results difficult to interpret. The spectra and apec fits are shown in Fig. 6, top right.

4.3 Region 3

Region 3 was well fitted by a single apec model with kT = 0.33 (0.24−0.67) keV with N_H = 0.08 (0.05−0.22) × 10²² cm⁻² giving a ${\chi }_{\nu }^2=1.06$ (see Table 2). As with Region 2, we also tried the nei model but the τ parameter was not well constrained so the model does not provide additional insight into the plasma conditions. The spectra and apec fits are shown in Fig. 6, bottom left.

4.4 Region 4

We fitted the spectra of Region 4 with a single apec model which provided a good fit to the data (${\chi }_{\nu }^2=0.98$) with a best-fit plasma temperature of 0.20 (0.18−0.21) keV and N_H < 0.04 × 10²² cm⁻² (see Table 2). No further constraints on the plasma conditions were provided using the nei model. The spectra and apec fits are shown in Fig. 6, bottom right.

Fig. 6 XMM-Newton spectral fits of our extended emission Regions 1 through 4. Each panel shows the EPIC spectrum with the most counts along with additive components of interest for clarity. The black data points and associated solid lines are theobservational spectra and best fit model, respectively. The magenta lines show the unresolved SNR, XRB, and SSS component (see Sect. 2.3.1 for descriptions), blue solid lines represent the M31 ISM component with the dashed blue lines indicating the contamination of smaller regions by larger extended emission regions, that is Region 5 components in Regions 2 and 3, and Region 7 component in Region 4 (see Fig. 8). Best fit parameters are given in Table 2. More detailed plots are presented in Fig. B.1.

Fig. 7 Same as Fig. 6 for Regions 5, 6, and 7, arranged from left to right. Best fit parameters are given in Table 2. More detailed plots are presented in Fig. B.2.

4.5 Region 5

The extended emission in Region 5 was well characterised by a single apec model (${\chi }_{\nu }^2=1.11$) with a best-fit plasma temperatureof 0.10 (0.07−0.15) keV and N_H = 0.18 (<0.34) × 10²² cm⁻² (see Table 2). The faintness of the extended emission in Region 5 over the background meant the nei model did not provide further information on the plasma conditions. The spectra and apec fits are shown in Fig. 7, left.

4.6 Region 6

A single apec model provided a good fit to the emission in Region 6 with ${\chi }_{\nu }^2=1.01$ for kT = 0.19 (0.18−0.22) keV and N_H < 0.39 × 10²² cm⁻² (see Table 2). We also applied the nei model which provided a constraint on the ionisation parameter with τ > 9.27 × 10¹¹ s cm⁻³ suggesting aplasma in collisional ionisation equilibrium. The spectra and nei fits are shown in Fig. 7, middle.

4.7 Region 7

Again, the extended emission in Region 7 was well accounted for by a single apec model (${\chi }_{\nu }^2=1.13$) with kT = 0.94 (0.86−1.99) keV and N_H = 0.20 (0.07−0.34) × 10²² cm⁻² (see Table 2). We also applied the nei model, though the fits did not yield any further information. The spectra and apec fits are shown in Fig. 7, right.

5 Discussion

5.1 Multi-wavelength view

The X-ray morphology is compared to H I, IR, Hα, and UV views of the northern disc in Fig. 8. The multi-wavelength view of the northern disc is dominated by the 10 kpc star-forming ring (Arp 1964; Habing et al. 1984). The prevalence of massive stars in the 10 kpc ring is clear from the Hα and FUV images (Fig. 8, bottom-left, bottom-right, and Fig. 9). The inner dust ring at ~ 5 kpc also hosts massive stellar clusters (MSCs) with significant UV flux there (Fig. 8, top-right, bottom-right).

The stellar cluster population of the northern disc of M31 has been catalogued using the PHAT data (Johnson et al. 2012, 2015). Ages and masses for many of these have been determined by Johnson et al. (2016). We used this catalogue to identify clusters with masses greater than 10³ M_⊙ with ages less than 100 Myr and plotted these over the NUV image, which is shown in Fig. 9. However, this study did not include the 10 kpc ring to the north and northwest of the M31 core so clusters in those regions are absent from Fig. 9. Nevertheless, there is an obvious correlation between the stellar cluster distribution and the extended X-ray emission, in particular along the 10 kpc ring to the northeast and east of the core, clear evidence of the link between the massive stellar population and the extended X-ray emission.

We found that the brightest extended X-ray emission in the northern disc follows closely the brightest regions in Hα and UV. Therefore, we conclude that the populations of massive stars in these regions play and important role in the production of the X-ray emission. These regions are discussed in more detail in Sect. 5.2. Fainter extended X-ray emission was detected in the low density regions between the 10 kpc and 5 kpc rings and between the 5 kpc ring and the core. The emission appears anti-correlated with the dust structures (Fig. 8, top-right) and does not appear to be associated with MSCs with intense UV flux as for Regions 1–4. These regions are discussed in detail in 5.3.

5.2 Emission from MSCs

5.2.1 Region 1

Region 1 spans what is probably the richest in massive stars region in the 10 kpc ring. Because of its extended nature along the ring itself, the extended X-ray emission likely results from several MSCs. This is evident from the multi-wavelength comparison shown in Fig. 10, top-left, with the extended X-ray emission clearly corresponding to enhancements in the Hα and NUV images. Knots of colder material, shown in the H I, 250 μm, and 24 μm images, appear distributed along the extent of Region 1, and could be sites of active star formation.

5.2.2 Region 2

The multi-wavelength view of Region 2 (Fig. 10, top-right) shows that X-ray emission likely results from several MSCs. Lewis et al. (2015) have shown that this region exhibited a high star formation rate in the last 25 Myr. The brightest region of the X-ray emission is centred on a large Hα shell measuring ~300 pc in diameter.The shell is bright on the eastern side though it appears broken in the west, possibly due to blowouts. Both the R-band and NUV images indicate an enhanced stellar density along the eastern edge of the shell, where the IR and H I emission is highest. The density gradient in these regions could explain the observed dichotomous morphology of the Hα shell. This multi-wavelength picture suggests that the X-ray emission in Region 2 could originate in a large superbubble. The size of the Hα shell is rather large for a superbubble, which typically have diameters of ~10² pc around young MSCs (e.g. Sasaki et al. 2011; Kavanagh et al. 2012; De Horta et al. 2014). It is possible that the Hα structure in Region 2 has been driven by more than one generation of MSCs, and is transitioning from a superbubble to a supergiant shell.

Fig. 8 Multi-wavelength comparison of hot ISM plasma in the northern disc of M31 with high mass stellar and cooler ISM components. In each case, the 0.3–0.7 keV X-ray image is given in blue and excised point sources have been filled. Top-left: comparison with H I (red). Top-right: comparison with the dust content, revealed in 250 μm (red) by Herschel SPIRE and 24 μm (green) by Spitzer MIPS. Bottom-left: Hα (red) from the LGGS, which traces the warm ISM; bottom-right: location of the massive stellar population, revealed in the NUV (red) and FUV (green) by GALEX. In all panels, the extended X-ray emission regions from our spectral analysis are shown to guide the eye, with labels given in top-right only.

5.2.3 Region 3

Region 3 was the brightest region of X-ray emission along the northern-most arc of the 10 kpc ring. The multi-wavelength view of Region 3 (Fig. 10, bottom-left) shows significant Hα emission towards the western side. The R-band and NUV images show MSCs are present, mostly concentrated in the centre of the region, with some extension to the southwest into the Hα structures. This region was found to have a high star formation rate over the last 50 Myr (Lewis et al. 2015). A slight east-west density gradient is also seen in the H I and IR images, which may explain the Hα morphology.

5.2.4 Region 4

The multi-wavelength view of Region 4 (see Fig. 10, bottom-right) shows no Hα structure associated with the X-ray enhancement, though the NUV image suggests that MSCs are present. Interestingly, the R-band image shows a slight deficit of stars compared to the surrounding regions. A density enhancement in H I, 250 μm, and 24 μm appears to be located in the foreground, absorbing the X-ray emission on the western side. This density enhancement may also explain the apparent deficit in stellar density seen in the R-band as resulting from foreground absorption.

5.2.5 Mechanical input from MSCs

As discussed earlier in this paper, the X-ray emission in Regions 1–4 most likely results from the mechanical input by massive stars and SNe in the MSCs of these regions. In this sub-section, we estimate the power injected into each region by the stellar populations using the PHAT photometric data with stellar evolution models and a population synthesis code. The results for each region are discussed below.

We constructed colour-magnitude diagrams for stars located within the regions using the F475W and F814W filters. To determine the properties of the stellar population we obtained stellar evolutionary tracks in the corresponding filters from the CMD 3.0 tool¹³, based on the PARSEC isochrones (Bressan et al. 2012; Chen et al. 2014, 2015; Tang et al. 2014). We selected the Chabrier lognormal IMF (Chabrier 2001) for M < 1 M_⊙ and the Salpeter IMF (Salpeter 1955) for M > 1 M_⊙, and assumed Galactic abundance. We extracted isochrones of various ages and evolutionary tracks at various masses and overplotted these on the PHAT photometry for the stellar populations in each region, shown in Fig. 11.

To assess the power input from the stellar population we use the Starburst 99 population synthesis code (Leitherer et al. 1999, 2010; Vázquez & Leitherer 2005). We estimated the total stellar mass of the initial population using the numbers of stars observed in the 8–12 M_⊙ region of the colour-magnitude diagrams and assumed an initial mass function (IMF) of the form used to calculate the PARSEC evolutionarytracks. We fed this total stellar mass into Starburst 99 and assumed an instantaneous burst of star formation to produce the stellar population. The input and output are detailed for each region below. The evolution of the power output due to stellar winds and supernovae was traced over time for each region and are shown in Fig. 12 along with the estimated age of the stellar populations. The stellar input properties for each region are summarised in Table 3.

Fig. 9 NUV image from Fig. 8 bottom-right, with stellar clusters from Johnson et al. (2016) shown by the red dots. The clusters have mass greater than 103 M⊙ and ages less than 100 Myr. The 0.4–1.25 keV X-ray contours from Fig. 2 are also shown in blue to highlight the correlation between the stellar clusters and X-ray emission. No clusters are plotted in the 10 kpc ring north of the core as this region was not studied in Johnson et al. (2016).

Region 1

Since Region 1 is quite elongated along the 10 kpc ring, we sub-divided the stellar population based on aggregations of the brightest stars found in the PHAT data. This resulted in four sub-populations which we label Regions 1a–1d. We constructed a colour-magnitude diagrams using the PHAT data (Fig. 11) and ran Starburst 99 simulations for each sub-region (Fig. 12). We found the most massive star was about 20–30 M_⊙ and ages ≲ 10 Myr for the MSCs. Using the Starburst 99 simulations, we estimated ~32−51 SNe have occurred in these MSCs, and they are currently injecting ~ 2.5−3.9 × 10³⁸ erg s⁻¹ into the region, with a total input energy over their lifetimes of ~ 6−8 × 10⁵² erg.

Region 2

For Region 2, the PHAT colour-magnitude diagram (Fig. 11) suggested the most massive star in the MSC is about 20 M_⊙ with an age of ~10 Myr for the MSC. The Starburst 99 simulations (Fig. 12) provided an estimate of ~ 15 past SNe in the region, with the MSC currently supplying ~ 5 × 10³⁷ erg s⁻¹ into the region, with integrated energy input of ~2 × 10⁵² erg.

Region 3

The PHAT colour-magnitude diagram (Fig. 11) suggested the most massive star in the MSC is ~ 30 M_⊙ with an MSC age of ~5−8 Myr. The Starburst 99 simulations (Fig. 12) showed that ~ 14−41 SNe may have already exploded in the region and the MSC is currently inputting ~ 2.0−4.0 × 10³⁸ erg s⁻¹ into the region, with integrated energy input of ~5−8 × 10⁵² erg.

Region 4

Finally, the PHAT photometry (Fig. 11) reveals the most massive star in the MSC to be ~ 20−30 M_⊙ and the age of the MSC to be ~8 Myr. Starburst99 calculated (Fig. 12) ~ 17 SNe have occurred and the MSC is currently pumping ~ 1.0 × 10³⁸ erg s⁻¹ into the region, with integrated energy input of ~3 × 10⁵² erg over the MSC lifetime.

5.2.6 Properties of the X-ray gas

We can estimate the physical properties of the hot gas in the extended emission region associated with MSCs using the morphology and the results of the X-ray spectral analysis (e.g. Sasaki et al. 2011; Kavanagh et al. 2012). We note that these estimates are somewhat coarse as we are limited in determining the dimensions of the emitting volumes because of the unknown depth of theemission regions in the line of sight. For Regions 2–4, which are approximately circular on the sky, we assume the volumes can be approximated as spherical. For Region 1, which is elongated along the 10 kpc ring, we assume the depth is the same as the projected width of the region (i.e. the volume is approximated as a highly elongated ellipsoid where the second and third principal axes are equal). While this is not ideal, the estimated plasma properties, though crude, can provide insight into the plasma conditions in the extended emission regions.

The gas density can be determined from the apec spectral model normalisation (K), which is defined as: $$K=\frac{1}{{10}^{14}\times 4\pi D_{\text{M} \text{31}}^2}{\displaystyle \int n_{\text{e}}}n \text{d}fV {\text{cm}}^{\text{−5}},$$(1)

where D_M 31 is the distance to M31, n_e is the electron density, n is the hydrogen density, and f and V are the filling parameter and volume of the X-ray emitting gas, respectively. We estimated the extent of the emission regions using the 5σ X-ray contours (see Fig. 10) and calculated a volume assuming a distance to M31 of 783 kpc (Conn et al. 2016). Assuming the metallicity of young clusters in M31 is approximately solar, the electron-to-hydrogen ratio n_e ∕n ≈ 1.21. Inserting the value K from Table 2 allowed us estimate the gas density for each region. The pressure (P_gas) in the gas can be estimated from the gas density and the temperature determined in the spectral fits (kT) as $$P_{\text{gas}}/k=2.31n f^{-1/2}T {\text{cm}}^{\text{−3}} \text{K}.$$(2)

Finally, the thermal energy contained in the X-ray emitting gas E can be estimated from $$E_{\text{gas}}=3/2 P_{\text{gas}}fV \text{erg}.$$(3)

The hot gas properties are summarised in Table 3.

The plasma properties determined for each of our regions, bearing in mind the caveat at the beginning of this subsection, are somewhat consistent with regions around MSCs and superbubbles of the well studied Magellanic Cloud (MC) populations. The gas densities in all regions are f^−1∕2 10⁻⁽²⁻³⁾ cm⁻³. Region 1 matches well with MC values (e.g. Sasaki et al. 2011; Kavanagh et al. 2012; Lopez et al. 2014; De Horta et al. 2014), though the density in Regions 2–4 is a factor of a few lower. This difference can most likely be explained by the difference in filling factor f between the MC objects and our M31 regions. Studies of MC objects are well resolved on pc scales so analyses can, for the most part, be focussed on individial MSCs and superbubbles. However, in the case of our regions, we cannot resolve X-ray emission from individual objects and we are observing the combined emission from several MSCs and superbubbles and the space between these meaning that the value of f is probably much lower than unity. A decrease in f causes a higher value of n with, e.g. f = 0.1 giving a factor of 3 increase in n. If this is indeed the case, gas densities in the M31 regions would be in close agreement with their MC counterparts. Similarly, the derived pressures are lower than expected. The pressures in regions of active star formation are thought to be P∕k = 10⁵⁻⁶ cm⁻³ K (Oey & García-Segura 2004, and references therein), and pressures in MSCs and stellar clusters should be higher. However, we find typical values P∕k of f^−1∕210⁴⁻⁵ cm⁻³ K in our regions. Again, if f is much lower than in the MC objects, as we suggest, this may explain the discrepancy. This should indeed be expected, since in a Milky Way like galaxy as M31, the level of star formation implies a supernova regulated ISM. It has been shown by 3D high resolution numerical simulations (de Avillez & Breitschwerdt 2004, 2005) that, e.g. in the solar neighbourhood, the volume filling factor f_h of hot gas with T > 10^5.5 K is f_h = 0.17 (increasing to 0.19 for a 3 μG field) for a Galactic supernova rate. This is a direct consequence of the high level of turbulence and a significant amount of gas being present in classically thermally unstable phases.

Other explanations for a decreased pressure in our regions are: (i) the possibility of significant blowouts into low density regions. The Hα-morphology of Region 2 exhibits a large ~300 pc diameter shell which appears to have blown out, possibly out of the 10 kpc ring. Such blowouts would result in a loss in pressure in the MSCs and superbubbles (see also Sect. 5.3.1); (ii) a natural cause for lower pressure is the deviation from CIE, which is often assumed in X-ray analyses, mostly for convenience. Also NEI fit models are only a crude representation of the real physics. The reason is because the ionisation structure has to be self-consistently calculated with the thermal one. The deviation from CIE is due to the fact that most X-ray plasmas are treated in the coronal approximation, where both collisional ionisation and X-ray emission extract energy from the plasma and are therefore cooling processes. So even if the plasma was in CIE initially, it will get out of CIE as it evolves in time (for details see Breitschwerdt & Schmutzler 1999; de Avillez & Breitschwerdt 2012). NEI will be especially enforced in regions of blowout, where adiabatic expansion of the plasma drastically reduces the kinetic energy of the electrons, while the recombination is delayed.

The determined thermal energies stored in the gas in our regions are slightly higher than expected with values of f^1∕2 10⁵²⁻⁵³ erg when compared to typical values in theMCs (10⁵¹⁻⁵² erg). This is not really surprising because, as mentioned above, each region contains several MSCs or superbubbles so the energy stored in a multitude of objects will naturally be higher than that in a single object as is usually studied in the MCs. The values of E_gas are also large when compared to the total energy input determined in the Starburst 99 simulations (E_SP,total in Table 3). This, once again, is most likely due to the real value of f which will act to reduce E_gas in the regions if f is much less than unity (see Eqs. (2) and (3)).

Fig. 10 Close-up, multi-wavelength views of extended emission regions asso- ciated with populations of massive stars with the wavebands indicated in the panes. The red, green, and blue contours correspond to 5, 10, and 20σ above the average background level determined from regions exterior to the northern disc. The individual regions are indicated onthe plot

Fig. 11 Colour-magnitude diagrams for the extended emission regions associated with MSCs. The data points are from PHAT photometry. Isochrones and evolutionary tracks derived from the PARSEC isochrone database are shown on the left and right, respectively.

Fig. 12 Wind and SN power input determined in Starburst 99 for the extended emission regions associated with MSCs. Wind powerand total energy evolution are shown by the green dashed lines, SNe are shown by the red dashed lines, and the combined input by the blue solid line. The red area show the age range of the stellar populations estimated from Fig. 11.

5.3 Large scale emission

5.3.1 Regions 5 and 6

The extended X-ray emission observed in Regions 5 and 6, at least in projection, appears to be located in the gap between the 5 kpc and 10 kpc rings in the northern disc (see Fig. 2) and is not obviously associated with any MSCs (see Fig. 8). The X-ray emission is soft in both regions (kT≲ 0.2 keV) and appearsto be close to or in CIE (see Table 2). In Region 5, the emission appears brightest near the 10 kpc ring MSCs in Region 3, suggesting that the MSCs play some role in this emission. The X-ray brightness drops off suddenly midway between the 5 kpc and 10 kpc rings in the south. It is unclear if this is a real reduction in X-ray brightness or is an absorption effect resulting from our viewing angle and the 5 kpc ring blocking our line of sight closer to its northern side. Certainly, such absorption effects can be invoked to explain why there is little X-ray emission observed in the gap between the rings near Region 1 where the dust 10 kpc ring is effectively masking the gap (Fig. 8, top-right).

Given the close correlation between the brightest emission in Region 5 and the MSCs in Region 3, and the nearby populations of massive stars in the 5 kpc ring for Region 6, it is possible that the MSCs in the ring are supplying the hot gas observed in the gaps. From Fig. 8, top-right and bottom-right, and Fig. 10, bottom right, we see that the MSCs in Region 3 are projected on the inner edge of the ring, and there is an ambient density gradient increasing away from the gap. Therefore, it is possible that the MSCs have heated the hot gas which has then blown out into the lower density region and into the gap. This may also be the case for the MSCs in the 5 kpc ring supplying the hot gas in Region 6, though it is possible that the location of these MSCs on the outer edge of the ring is caused by projection. An alternative explanation for the X-ray emission observed in Region 6 is that the hot gas is a relic of an older massive stellar population, which no longer hosts the massive stars responsible for the heating. However, the cooling time for the hot plasma is ~ 10⁶⁻⁷ yr, so additional heating by the MSCs in the 5 kpc ring may still be required in this case.

5.3.2 Region 7

The extended emission in Region 7 extends from the M31 core out to the 5 kpc ring (excluding the emission in Region 4). This region encompasses both gaps between dust structures inside the 5 kpc ring (Fig. 8, top-right) and also some populations of massive stars (Fig. 8, bottom-right), though comparatively fewer than in the 10 kpc ring and still anti-correlated with the brightest X-ray emission. The X-ray emission is harder than in other extended regions, including those associated with MSCs and superbubbles in the 10 kpc ring, with kT = 0.94 (0.86 − 0.99) keV compared to kT ≲ 0.7 keV. It is unclear as to why this may be, though it is obvious that Region 7 contains a much more diverse multi-wavelength morphology and stellar population than Regions 6 and 7, so some additional heating mechanism or a longer cooling time due to lower density may be responsible. In addition, it could also be that there is some outflow component originating in Region 4, as seen in NGC 5461 in M 101 for example (Sun et al. 2012), projected against Region 7. The observed emission is likely an amalgam of several different emission regions and mechanisms. The low count statistics makes it difficult for us to sub-divide the region and perform a more meaningful spectral analysis using smaller spatial scales.

5.4 Comparison with other spiral galaxies

Early observations of nearby galaxies in X-rays with the Einstein Observatory and ROSAT revealed soft diffuse X-rays fromnearby spiral galaxy (Fabbiano et al. 1992; Snowden & Pietsch 1995; Ehle et al. 1995; Cui et al. 1996; Read et al. 1997). Newer observations with XMM-Newton and Chandra have shown that the diffuse X-ray emission in galaxy discs is found along the spiral arms and is correlated with regions with enhanced Hα, UV, or IR emission (Kuntz et al. 2003; Tyler et al. 2004). The multi-wavelength correlation indicates that soft diffuse X-ray emission traces the regions in the disc with active star formation. The morphology and the spectrum of the diffuse emission was studied in detail in, e.g. M 101 using Chandra and XMM-Newton data by Kuntz et al. (2003) and Warwick et al. (2007), respectively.The spectrum of the diffuse X-ray emission in the disc of spiral galaxies can be well reproduced by a thermal model with kT = 0.20 keV, while in some cases, an additional component with a higher temperature (kT= 0.6−0.8 keV) is required. The soft component can be explained as emission from the thermal hot phase of the ISM, while the hotter component can have contribution of late-type stars or binaries (see also Sect. 2.3.2). We studied the nearby pair of galaxies NGC 1512 and NGC 1510 with XMM-Newton. While the photon statistics were not high enough to detect diffuse emission in the outer parts of the disc, we detected significant diffuse emission in the central region of the main galaxy NGC 1512 with kT = 0.66 keV (Ducci et al. 2014).

In the northern disc of M31, significant diffuse X-ray emission is detected along the major dust ring, with some additional diffuse X-ray emission filling the region inside the dust ring. Regions with bright diffuse emission seem to correlate with regions with active star formation as seen in the distribution of young massive stars in UV, H II regions in Hα, or dust in IR. The temperature of the hot plasma varies from kT = 0.1–0.3 keV in the ring up to kT = 0.6 keV in a superbubble. The region inside the major dust ring and along the inner rings appears to be filled with thermal plasma at a higher temperature (kT = 0.94 keV, Region 7) even though there are no obvious bright star-forming regions. Together with the disturbed appearance of the ring which consists of narrower fainter sub-structures, the higher temperature might indicate that the ISM was not only heated by massive stars but also by some disturbances in the disc. Overall, with a few notable exceptions, the morphology and spectra of the hot gas in the northern disc of M31 is similar to that of other galaxy discs.

6 Conclusion

We used new deep and archival XMM-Newton observations of the northern disc of M31 to trace the hot ISM plasma in unprecedented detail and to characterise the physical properties of the X-ray emitting gas. We also used multi-wavelength surveys to put the X-ray emission in context of the multi-phase ISM.

The brightest extended emission regions were found to be correlated with regions hosting populations of massive stars, notably in the 5 kpc and 10 kpc star-forming rings. The plasma temperatures in these regions were ~ 0.2 up to ~0.6 keV, which are consistent with an origin in populations of massive stars and superbubbles. The derived X-ray luminosities, densities, and pressures for the gas in each of these regions are all consistent with typical values in the literature for the well-studied Magellanic Cloud populations if the filling factor of the emission is much lower than unity. We argue that this is indeed the case for the regions studied since we don’t have the spatial resolution to disentangle individual MSCs and superbubbles.

We found large extended emission filling low density gaps in the dust morphology of the northern disc, most notably between the 5 and 10 kpc star-forming rings. We suggested that the hot gas in the low density regions most likely originated in the MSCs in the rings and expelled, moving down the density gradient into the gaps.

Our deep observations of the northern disc have demonstrated that the hot ISM plasma in a galaxy like our own can be studied on relatively small spatial scales without the issues such as foreground absorption, distance uncertainties, and source confusion that affect Galactic studies. The southern disc of M31 offers another promising target for such a deep study. The 10 kpc ring in the southern disc hosts the large stellar association NGC 206, and is located where the ring was most likely disturbed by an encounter with the satellite galaxy M 32. Therefore, deep observations of the southern disc in X-rays would allow the study of, not only the hot ISM plasma associated with massive stellar clusters, but also the environment of the M31, M 32 collision region.

Acknowledgements

We wish to thank the anonymous referee for their comments and suggestions which improved the paper. This research is funded by the German Bundesministerium für Wirtschaft und Technologie and the Deutsches Zentrum für Luft- und Raumfahrt (BMWi/DLR) through the grant FKZ 50 OR 1510. M.S. acknowledges support by the Deutsche Forschungsgemeinschaft through the Heisenberg professor grants SA 2131/5-1 and 12-1. This research made use of Montage. It is funded by the National Science Foundation under Grant Number ACI-1440620, and was previously funded by the National Aeronautics and Space Administration’s Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology.

Appendix A Table of observations

Appendix B Spectral fits

Fig. B.1 XMM-Newton spectral fits of our extended emission Regions 1 through 4, arranged from top to bottom. Left panels: fits to all spectra. Right panels: EPIC spectra with the most counts along with additive components of the best fit model for clarity. The red data points and line are the FWC data and best fit. The black data points and associated lines are theobservational spectra, best fit model (solid), combined astrophysical components (dash-dot), and combined particle-induced components (FWC plus residual SP, dashed). The orange lines represent the AXB components, the magenta line shows the unresolved SNR, XRB, and SSS component (see Sect. 2.3.1 for descriptions), blue solid lines represent the M31 ISM component with the dashed blue lines indicating the contamination of smaller regions by larger extended emission regions (e.g. Region 5 component in Region 2). Best fit parameters are given in Table 2.

Fig. B.2 Same as Fig. B.1 for Regions 5, 6, and 7, arranged from top to bottom. Best fit parameters are given in Table 2.

References

All Tables

All Figures

Fig. 1 XMM-Newton mosaic images of M31 with 0.4–0.7 keV shown in top-left and 0.7–1.25 keV in top-right. The white dashed circles mark the fields of our northern disc LP survey, which are revealed in more detail in the zoomed images in bottom-left and bottom-right. The XMM-Newton coverage in relation to the stellar population is shown in Fig. 1 of Stiele et al. (2011). The colours are in linear scale as shown by the colour bars in each image. The units of the images are cts s−1 deg−2.

In the text

Fig. 2 Left: XMM-Newton 0.4–1.25 keV image of the northern disc of M31 with contours. The image has been adaptively smoothed and masked point sources filled as described in the text. Yellow, cyan, blue, and red contours correspond to 5, 10, 20, and 50σ above the average background level determined from regions exterior to the northern disc. Regions of interest discussed in the text are marked. Right: XMM-Newton 0.4–1.25 keV image of the northern disc with the spectral extraction regions defined from the contours in left shown in cyan. Background regions are shown in yellow. The deep EPIC northern disc fields are shown by the green circles.

In the text

Fig. 3 XMM-Newton spectral fits of the background regions. Top-left: simultaneous spectral fits to the background region NE, and corresponding FWC spectra extracted from EPIC-MOS2 (green and blue) and EPIC-pn (black and red). The dashed lines show the combined AXB contribution while the dotted lines show the combined particle-induced contribution to the extended emission spectra. For clarity, only spectra from ObsID 0763120301 are shown. Top-right: same as top-left but for background region NW. Bottom-left: same as top-left but for background region SE. In this case the spectra from ObsID 0763120101 are shown. Bottom-right: simultaneous fits to all background spectra with each background region source and FWC spectrum shown in a different colour. The combined AXB and FWC contributions are also shown for each. See text for details.

In the text

Fig. 4 X-ray flux distributions of M31 sources detected in the southern (top row) and northern (bottom row) deep fields in the northern disc. Left panels: X-ray flux functions, along with power law and cut-off power law fits. Right panels: combined flux distribution expected from sources in flux bins, normalised for the deep field areas. The unresolved flux per unit area is determined by integrating over this distribution and subtracting the flux of detected sources.

In the text

	Fig. 5 Combined faint source spectrum, fitted with an absorbed power law model with best fit NH = 0.96 (0.74−1.2) × 1021 cm−2 and Γ = 2.37 (2.22−2.53) giving a reduced χ2 = 1.07 for 662 degrees of freedom.
In the text

Fig. 6 XMM-Newton spectral fits of our extended emission Regions 1 through 4. Each panel shows the EPIC spectrum with the most counts along with additive components of interest for clarity. The black data points and associated solid lines are theobservational spectra and best fit model, respectively. The magenta lines show the unresolved SNR, XRB, and SSS component (see Sect. 2.3.1 for descriptions), blue solid lines represent the M31 ISM component with the dashed blue lines indicating the contamination of smaller regions by larger extended emission regions, that is Region 5 components in Regions 2 and 3, and Region 7 component in Region 4 (see Fig. 8). Best fit parameters are given in Table 2. More detailed plots are presented in Fig. B.1.

In the text

	Fig. 7 Same as Fig. 6 for Regions 5, 6, and 7, arranged from left to right. Best fit parameters are given in Table 2. More detailed plots are presented in Fig. B.2.
In the text

Fig. 8 Multi-wavelength comparison of hot ISM plasma in the northern disc of M31 with high mass stellar and cooler ISM components. In each case, the 0.3–0.7 keV X-ray image is given in blue and excised point sources have been filled. Top-left: comparison with H I (red). Top-right: comparison with the dust content, revealed in 250 μm (red) by Herschel SPIRE and 24 μm (green) by Spitzer MIPS. Bottom-left: Hα (red) from the LGGS, which traces the warm ISM; bottom-right: location of the massive stellar population, revealed in the NUV (red) and FUV (green) by GALEX. In all panels, the extended X-ray emission regions from our spectral analysis are shown to guide the eye, with labels given in top-right only.

In the text

Fig. 9 NUV image from Fig. 8 bottom-right, with stellar clusters from Johnson et al. (2016) shown by the red dots. The clusters have mass greater than 103 M⊙ and ages less than 100 Myr. The 0.4–1.25 keV X-ray contours from Fig. 2 are also shown in blue to highlight the correlation between the stellar clusters and X-ray emission. No clusters are plotted in the 10 kpc ring north of the core as this region was not studied in Johnson et al. (2016).

In the text

	Fig. 10 Close-up, multi-wavelength views of extended emission regions asso- ciated with populations of massive stars with the wavebands indicated in the panes. The red, green, and blue contours correspond to 5, 10, and 20σ above the average background level determined from regions exterior to the northern disc. The individual regions are indicated onthe plot
In the text

	Fig. 11 Colour-magnitude diagrams for the extended emission regions associated with MSCs. The data points are from PHAT photometry. Isochrones and evolutionary tracks derived from the PARSEC isochrone database are shown on the left and right, respectively.
In the text

	Fig. 12 Wind and SN power input determined in Starburst 99 for the extended emission regions associated with MSCs. Wind powerand total energy evolution are shown by the green dashed lines, SNe are shown by the red dashed lines, and the combined input by the blue solid line. The red area show the age range of the stellar populations estimated from Fig. 11.
In the text

Fig. B.1 XMM-Newton spectral fits of our extended emission Regions 1 through 4, arranged from top to bottom. Left panels: fits to all spectra. Right panels: EPIC spectra with the most counts along with additive components of the best fit model for clarity. The red data points and line are the FWC data and best fit. The black data points and associated lines are theobservational spectra, best fit model (solid), combined astrophysical components (dash-dot), and combined particle-induced components (FWC plus residual SP, dashed). The orange lines represent the AXB components, the magenta line shows the unresolved SNR, XRB, and SSS component (see Sect. 2.3.1 for descriptions), blue solid lines represent the M31 ISM component with the dashed blue lines indicating the contamination of smaller regions by larger extended emission regions (e.g. Region 5 component in Region 2). Best fit parameters are given in Table 2.

In the text

	Fig. B.2 Same as Fig. B.1 for Regions 5, 6, and 7, arranged from top to bottom. Best fit parameters are given in Table 2.
In the text