QUIJOTE Scientific Results – XVII. Studying the Anomalous Microwave Emission in the Andromeda Galaxy with QUIJOTE-MFI

M Fernández-Torreiro, R.T Génova-Santos, J.A Rubiño-Martín, C.H López-Caraballo, M.W Peel, C Arce-Tord, R Rebolo, E Artal, M Ashdown, R.B Barreiro, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

M Fernández-Torreiro, R.T Génova-Santos, J.A Rubiño-Martín, C.H López-Caraballo, M.W Peel, et al.. QUIJOTE Scientific Results – XVII. Studying the Anomalous Microwave Emission in the Andromeda Galaxy with QUIJOTE-MFI. Monthly Notices of the Royal Astronomical Society, 2024, 527 (4), pp.11945-11961. 10.1093/mnras/stad3145. hal-04103424

HAL Id: hal-04103424

https://hal.science/hal-04103424

Submitted on 6 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# QUIJOTE Scientific Results – XVII. Studying the anomalous microwave emission in the Andromeda Galaxy with QUIJOTE-MFI

M. Fernández-Torreiro<sup>1,1,2*</sup> R. T. Génova-Santos<sup>1,1,2*</sup> J. A. Rubiño-Martín<sup>1,2</sup>

C. H. López-Caraballo<sup>1,2</sup> M. W. Peel<sup>1,2,3</sup> C. Arce-Tord<sup>4,5</sup> R. Rebolo<sup>1,2,6</sup> E. Artal<sup>7</sup>

M. Ashdown<sup>8,9</sup> R. B. Barreiro<sup>10</sup> F. J. Casas<sup>10</sup> E. de la Hoz<sup>10,11,12</sup> F. Guidi<sup>13</sup> D. Herranz<sup>10</sup>

R. Hoyland<sup>1,2</sup> A. Lasenby<sup>8,9</sup> E. Martínez-Gonzalez<sup>10</sup> L. Piccirillo,<sup>14</sup> F. Poidevin<sup>1,2</sup>

B. Ruiz-Granados<sup>1,2,15</sup> D. Tramonte,<sup>16,1,2</sup> F. Vansyngel<sup>1,2</sup> P. Vielva<sup>10</sup> and R. A. Watson<sup>14</sup>

*Affiliations are listed at the end of the paper*

Accepted 2023 October 10. Received 2023 October 10; in original form 2023 May 12

## ABSTRACT

The Andromeda Galaxy (M31) is the Local Group galaxy that is most similar to the Milky Way (MW). The similarities between the two galaxies make M31 useful for studying integrated properties common to spiral galaxies. We use the data from the recent QUIJOTE-MFI Wide Survey, together with new raster observations focused on M31, to study its integrated emission. The addition of raster data improves the sensitivity of QUIJOTE-MFI maps by almost a factor 3. Our main interest is to confirm if anomalous microwave emission (AME) is present in M31, as previous studies have suggested. To do so, we built the integrated spectral energy distribution of M31 between 0.408 and 3000 GHz. We then performed a component separation analysis taking into account synchrotron, free-free, AME, and thermal dust components. AME in M31 is modelled as a log-normal distribution with maximum amplitude, $A_{\text{AME}}$, equal to $1.03 \pm 0.32$ Jy. It peaks at $v_{\text{AME}} = 17.2 \pm 3.2$ GHz with a width of $W_{\text{AME}} = 0.58 \pm 0.16$. Both the Akaike and Bayesian information criteria find the model without AME to be less than 1 per cent as probable as the one taking AME into consideration. We find that the AME emissivity per 100 µm intensity in M31 is $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 9.6 \pm 3.1$ µK MJy$^{-1}$ sr, similar to that of the MW. We also provide the first upper limits for the AME polarization fraction in an extragalactic object. M31 remains the only galaxy where an AME measurement has been made of its integrated spectrum.

**Key words:** radiation mechanisms: general – ISM: general – galaxies: ISM – Local Group – diffuse radiation – radio continuum: galaxies.

## 1 INTRODUCTION

The Andromeda Galaxy or Messier 31 (M31) is the largest galaxy in the Local Group (Goodwin, Gribbin & Hendry 1998) and the most similar to the Milky Way (MW), both being spiral galaxies. It has been observed and studied in detail throughout modern astronomical history owing to its large angular size,<fnref n="1" /> over 5 deg². These studies span the whole frequency domain: a review of high-resolution radio measurements is available in Berkhuijsen, Beck & Hoernes (2003). Fatigoni et al. (2021) also studied M31 at high-angular resolution, using data taken in the C-band of the Sardinia Radio Telescope (hereafter, SRT). Together with the previously mentioned ancillary data, thermal and non-thermal emission components

could be spatially disentangled, and the star formation rate was estimated.

Anomalous microwave emission (hereafter, AME) is a prominent mechanism of Galactic emission in the frequency range 10–60 GHz. Known for more than 25 years now (Kogut et al. 1996; Leitch et al. 1997), the physical process responsible for AME is not clear yet. The dominant hypothesis states that spinning dust emission is responsible for AME (Draine & Lazarian 1998a, b), owing partly to the strong spatial correlation between AME and thermal dust (300–3000 GHz) emission (e.g. Fernández-Torreiro et al. 2023; Poidevin et al. 2023). However, the exact carrier responsible for AME is still a mystery. Neither are its polarization properties well understood yet, although observations both on compact regions at angular scales of 1° or below (e.g. Dickinson, Peel & Vidal 2011; López-Caraballo et al. 2011; Génova-Santos et al. 2015, 2017) and on large angular scales (Macellari et al. 2011; Herman et al. 2023) indicate that its polarization fraction should be below 5 per cent. See Dickinson et al. (2018) for a detailed review of AME.

<footnote>*E-mail: mateo.fernandez@iac.es (MFT); rgs@iac.es (RTGS)</footnote>

<footnote n="1">¹M31’s isophotal major radius is 91.5 arcmin (de Vaucouleurs et al. 1991) with an axial ratio of ≈0.7. The Herschel Exploitation of Local Galaxy Andromeda (HELGA, Fritz et al. 2012) program went beyond M31 itself and also studied its surroundings, producing maps of ~5.5° × 2.5°.</footnote>
# QUIJOTE Scientific Results – XVII. Studying the anomalous microwave emission in the Andromeda Galaxy with QUIJOTE-MFI

M. Fernández-Torreiro<sup>1,1,2*</sup> R. T. Génova-Santos<sup>1,1,2*</sup> J. A. Rubiño-Martín<sup>1,2</sup>

C. H. López-Caraballo<sup>1,2</sup> M. W. Peel<sup>1,2,3</sup> C. Arce-Tord<sup>4,5</sup> R. Rebolo<sup>1,2,6</sup> E. Artal<sup>7</sup>

M. Ashdown<sup>8,9</sup> R. B. Barreiro<sup>10</sup> F. J. Casas<sup>10</sup> E. de la Hoz<sup>10,11,12</sup> F. Guidi<sup>13</sup> D. Herranz<sup>10</sup>

R. Hoyland<sup>1,2</sup> A. Lasenby<sup>8,9</sup> E. Martínez-Gonzalez<sup>10</sup> L. Piccirillo,<sup>14</sup> F. Poidevin<sup>1,2</sup>

B. Ruiz-Granados<sup>1,2,15</sup> D. Tramonte,<sup>16,1,2</sup> F. Vansyngel<sup>1,2</sup> P. Vielva<sup>10</sup> and R. A. Watson<sup>14</sup>

*Affiliations are listed at the end of the paper*

Accepted 2023 October 10. Received 2023 October 10; in original form 2023 May 12

## ABSTRACT

The Andromeda Galaxy (M31) is the Local Group galaxy that is most similar to the Milky Way (MW). The similarities between the two galaxies make M31 useful for studying integrated properties common to spiral galaxies. We use the data from the recent QUIJOTE-MFI Wide Survey, together with new raster observations focused on M31, to study its integrated emission. The addition of raster data improves the sensitivity of QUIJOTE-MFI maps by almost a factor 3. Our main interest is to confirm if anomalous microwave emission (AME) is present in M31, as previous studies have suggested. To do so, we built the integrated spectral energy distribution of M31 between 0.408 and 3000 GHz. We then performed a component separation analysis taking into account synchrotron, free-free, AME, and thermal dust components. AME in M31 is modelled as a log-normal distribution with maximum amplitude, $A_{\text{AME}}$, equal to $1.03 \pm 0.32$ Jy. It peaks at $v_{\text{AME}} = 17.2 \pm 3.2$ GHz with a width of $W_{\text{AME}} = 0.58 \pm 0.16$. Both the Akaike and Bayesian information criteria find the model without AME to be less than 1 per cent as probable as the one taking AME into consideration. We find that the AME emissivity per 100 µm intensity in M31 is $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 9.6 \pm 3.1$ µK MJy$^{-1}$ sr, similar to that of the MW. We also provide the first upper limits for the AME polarization fraction in an extragalactic object. M31 remains the only galaxy where an AME measurement has been made of its integrated spectrum.

**Key words:** radiation mechanisms: general – ISM: general – galaxies: ISM – Local Group – diffuse radiation – radio continuum: galaxies.

## 1 INTRODUCTION

The Andromeda Galaxy or Messier 31 (M31) is the largest galaxy in the Local Group (Goodwin, Gribbin & Hendry 1998) and the most similar to the Milky Way (MW), both being spiral galaxies. It has been observed and studied in detail throughout modern astronomical history owing to its large angular size,<fnref n="1" /> over 5 deg². These studies span the whole frequency domain: a review of high-resolution radio measurements is available in Berkhuijsen, Beck & Hoernes (2003). Fatigoni et al. (2021) also studied M31 at high-angular resolution, using data taken in the C-band of the Sardinia Radio Telescope (hereafter, SRT). Together with the previously mentioned ancillary data, thermal and non-thermal emission components

could be spatially disentangled, and the star formation rate was estimated.

Anomalous microwave emission (hereafter, AME) is a prominent mechanism of Galactic emission in the frequency range 10–60 GHz. Known for more than 25 years now (Kogut et al. 1996; Leitch et al. 1997), the physical process responsible for AME is not clear yet. The dominant hypothesis states that spinning dust emission is responsible for AME (Draine & Lazarian 1998a, b), owing partly to the strong spatial correlation between AME and thermal dust (300–3000 GHz) emission (e.g. Fernández-Torreiro et al. 2023; Poidevin et al. 2023). However, the exact carrier responsible for AME is still a mystery. Neither are its polarization properties well understood yet, although observations both on compact regions at angular scales of 1° or below (e.g. Dickinson, Peel & Vidal 2011; López-Caraballo et al. 2011; Génova-Santos et al. 2015, 2017) and on large angular scales (Macellari et al. 2011; Herman et al. 2023) indicate that its polarization fraction should be below 5 per cent. See Dickinson et al. (2018) for a detailed review of AME.

<footnote>*E-mail: mateo.fernandez@iac.es (MFT); rgs@iac.es (RTGS)</footnote>

<footnote n="1">¹M31’s isophotal major radius is 91.5 arcmin (de Vaucouleurs et al. 1991) with an axial ratio of ≈0.7. The Herschel Exploitation of Local Galaxy Andromeda (HELGA, Fritz et al. 2012) program went beyond M31 itself and also studied its surroundings, producing maps of ~5.5° × 2.5°.</footnote>
While there have been many detections of AME in our Galaxy (e.g. Watson et al. 2005; Planck Collaboration et al. 2014a; Poidevin et al. 2019; Cepeda-Arroita et al. 2021; Fernández-Torreiro et al. 2023) extragalactic detections of AME are scarce (Murphy et al. 2010, 2018). Peel et al. (2011) set upper limits on the AME from three different bright galaxies using the WMAP 7 yr Point Source Catalogue and the Planck Early Release Compact Source Catalogue (ERCSC). Using observations from these two satellites, Planck Collaboration et al. (2015b) claimed a 2.3σ significance detection of AME in M31. Tibbs et al. (2018) also analysed Planck data for M33 but found no evidence of an AME component on its integrated spectrum. Bianchi et al. (2022) studied data from K-band detectors at SRT for four nearby spiral galaxies but was unable to detect AME in any of them, which instead places constraints on its upper levels. However, AME has been detected in resolved extragalactic regions, such as a star-forming region in NGC 6946 (Murphy et al. 2010; Hensley, Murphy & Staguhn 2015) and a compact radio source associated with NGC 4725 (Murphy et al. 2018). It appears that AME is much easier to detect in individual, well isolated regions, where the background emission can be easily subtracted. Generally, these regions also show high emission from thermal dust and high star formation ratios, together with low synchrotron signal. It is difficult to find a galaxy whose integrated spectrum shows this combination of properties, as normally many strong components contribute to the spectrum.

Based on new data at 6.6 GHz from the SRT, which are key to determining the level of free-free emission, Battistelli et al. (2019) (hereafter, B19) claimed a high-significance detection of AME in M31. These SRT data were smoothed to 1° resolution and studied together with data from WMAP, Planck and Herschel (Fritz et al. 2012) satellites. Few data sets between 2 and 20 GHz covering M31 are currently available. This causes synchrotron, free-free and AME models to be highly degenerate unless strong priors are placed on their parameters. However, the combination of SRT data with lower frequency (below 2 GHz) surveys permitted a good determination of the synchrotron and free-free levels. Once synchrotron and free-free have been properly defined, a good determination of AME is straightforward. The significance of the AME detection increased from 2.3σ in Planck Collaboration et al. (2015b) to over 8σ in B19 after including these SRT data.

Harper et al. (2023) performed an independent analysis of M31 after adding data from the C-BASS experiment at 4.76 GHz. This data addition should improve the disentanglement of synchrotron and free-free, similar to what happened with SRT data. Using the same apertures as Planck Collaboration et al. (2015b), they claimed a 3.0σ detection of AME, although with much lower amplitude than that expected from both Planck Collaboration et al. (2015b) and B19. Planck Collaboration et al. (2015b) and Harper et al. (2023) obtained $S_{\text{AME}}^{30\text{GHz}} = 0.7 \pm 0.3$ and $0.27 \pm 0.09$ Jy, respectively, and B19 obtained $S_{\text{AME}}^{25\text{GHz}} = 1.45 \pm 0.15$ Jy while using a smaller aperture. As a consequence of their low AME amplitude, Harper et al. (2023) also obtained an extremely low estimate of the AME emmissivity, a factor 20 lower than those expected from measurements on our Galaxy (Planck Collaboration et al. 2015a; Harper et al. 2022; Fernández-Torreiro et al. 2023).

In this paper, we present the data taken on M31 with the Q-U-I JOint Tenerife Experiment Multi Frequency Instrument (QUIJOTE-MFI) at 11 and 13 GHz between 2012 and 2018 from the combination between its Wide Survey (Rubiño-Martín et al. 2023) and focused raster observations. We use these QUIJOTE-MFI data to assess the presence of AME in M31 following an independent procedure to that of B19 and Harper et al. (2023). The paper is organized as follows.

Section 2 describes the data used to build the M31 integrated spectral energy distribution (SED) and how they were processed. Section 3 describes the components assumed in the fitting of the SED of M31 and the fitting procedure itself. Section 4 presents the main results of this study, while Section 5 comments on the changes for these results with different assumptions. Finally, we present the main conclusions of this work in Section 6.

## 2 INPUT DATA

A summary of the maps used is presented in Table 1. All of them are used in HEALPix<fnref n="2" /> (Górski et al. 2005; Zonca et al. 2019) pixelization $N_{\text{side}} = 512$ and smoothed to a common resolution of 1°.

### 2.1 QUIJOTE-MFI data

The QUIJOTE CMB experiment (Rubiño-Martín et al. 2010) operates at the Teide Observatory (OT) of the Instituto de Astrofísica de Canarias (IAC), located at latitude $28^\circ 18'04''$ North and longitude $16^\circ 30'38''$ West. The first instrument installed on QUIJOTE was the Multi Frequency Instrument (MFI). This instrument observed the sky, in intensity and polarization, at 11, 13, 17, and 19 GHz, mostly in the so-called 'nominal mode' configuration, with continuous observations at constant elevation while spinning in azimuth. These observations constitute the Wide Survey (hereafter, WS), which is described in detail in Rubiño-Martín et al. (2023). In addition to the WS data set, we use new data obtained with raster mode observations specifically focused on M31, intended to improve the sensitivity in this area. All the analyses and results presented in this paper have been derived using 11 and 13 GHz maps. High noise levels in the 17 and 19 GHz intensity maps (see Table 2), mostly dictated by atmospheric 1/f noise, prevent the detection of emission from M31 at these frequencies so they are discarded for this analysis. The 11 and 13 GHz maps are produced combining the WS data with the data obtained in raster scan mode around the coordinates of M31. A single run of the QUIJOTE Map-Making code (Guidi et al. 2021) is required in order to optimize the recovering of the large angular-scale signal. The same strategy has been applied before (Guidi et al. 2023; Tramonte et al. 2023) and will be used in upcoming (González-González et al., in preparation; Ruiz-Granados et al., in preparation) publications using QUIJOTE-MFI data.

The data in raster-scan mode were taken between May and December 2016 for a total of 931 rasters accounting for 539.1 h of observations (see Table 3). These raster-scans are constant-elevation azimuth scans of width $\approx 12^\circ/\cos(\text{EL})$ and duration $\approx 35$ min, performed on different local coordinates while tracking the field, with elevations ranging between $32^\circ$ and $75^\circ$. After projection onto the sky plane, each of these observations results in a map of $\approx 12^\circ$ by $\approx 12^\circ$. The combination of all these observations results in a sky coverage as shown in the bottom-right-hand panel of Fig. 1, where the footprint of the raster scans is evident. This is one of the fields with highest integration time with the MFI. As a result, these are the maps with the best sensitivity of all of those obtained with the MFI to date, both in intensity and polarization. Final map sensitivities, computed through a jack-knife analysis on the Half Mission Difference Maps (HMDM; see Rubiño-Martín et al. 2023) are shown in Table 2: we can clearly see the improvement in sensitivity owing to the addition of raster data for most of the horns. In intensity (which are the data

<footnote n="2">²https://healpix.sourceforge.io</footnote>
While there have been many detections of AME in our Galaxy (e.g. Watson et al. 2005; Planck Collaboration et al. 2014a; Poidevin et al. 2019; Cepeda-Arroita et al. 2021; Fernández-Torreiro et al. 2023) extragalactic detections of AME are scarce (Murphy et al. 2010, 2018). Peel et al. (2011) set upper limits on the AME from three different bright galaxies using the WMAP 7 yr Point Source Catalogue and the Planck Early Release Compact Source Catalogue (ERCSC). Using observations from these two satellites, Planck Collaboration et al. (2015b) claimed a 2.3σ significance detection of AME in M31. Tibbs et al. (2018) also analysed Planck data for M33 but found no evidence of an AME component on its integrated spectrum. Bianchi et al. (2022) studied data from K-band detectors at SRT for four nearby spiral galaxies but was unable to detect AME in any of them, which instead places constraints on its upper levels. However, AME has been detected in resolved extragalactic regions, such as a star-forming region in NGC 6946 (Murphy et al. 2010; Hensley, Murphy & Staguhn 2015) and a compact radio source associated with NGC 4725 (Murphy et al. 2018). It appears that AME is much easier to detect in individual, well isolated regions, where the background emission can be easily subtracted. Generally, these regions also show high emission from thermal dust and high star formation ratios, together with low synchrotron signal. It is difficult to find a galaxy whose integrated spectrum shows this combination of properties, as normally many strong components contribute to the spectrum.

Based on new data at 6.6 GHz from the SRT, which are key to determining the level of free-free emission, Battistelli et al. (2019) (hereafter, B19) claimed a high-significance detection of AME in M31. These SRT data were smoothed to 1° resolution and studied together with data from WMAP, Planck and Herschel (Fritz et al. 2012) satellites. Few data sets between 2 and 20 GHz covering M31 are currently available. This causes synchrotron, free-free and AME models to be highly degenerate unless strong priors are placed on their parameters. However, the combination of SRT data with lower frequency (below 2 GHz) surveys permitted a good determination of the synchrotron and free-free levels. Once synchrotron and free-free have been properly defined, a good determination of AME is straightforward. The significance of the AME detection increased from 2.3σ in Planck Collaboration et al. (2015b) to over 8σ in B19 after including these SRT data.

Harper et al. (2023) performed an independent analysis of M31 after adding data from the C-BASS experiment at 4.76 GHz. This data addition should improve the disentanglement of synchrotron and free-free, similar to what happened with SRT data. Using the same apertures as Planck Collaboration et al. (2015b), they claimed a 3.0σ detection of AME, although with much lower amplitude than that expected from both Planck Collaboration et al. (2015b) and B19. Planck Collaboration et al. (2015b) and Harper et al. (2023) obtained $S_{\text{AME}}^{30\text{GHz}} = 0.7 \pm 0.3$ and $0.27 \pm 0.09$ Jy, respectively, and B19 obtained $S_{\text{AME}}^{25\text{GHz}} = 1.45 \pm 0.15$ Jy while using a smaller aperture. As a consequence of their low AME amplitude, Harper et al. (2023) also obtained an extremely low estimate of the AME emmissivity, a factor 20 lower than those expected from measurements on our Galaxy (Planck Collaboration et al. 2015a; Harper et al. 2022; Fernández-Torreiro et al. 2023).

In this paper, we present the data taken on M31 with the Q-U-I JOint Tenerife Experiment Multi Frequency Instrument (QUIJOTE-MFI) at 11 and 13 GHz between 2012 and 2018 from the combination between its Wide Survey (Rubiño-Martín et al. 2023) and focused raster observations. We use these QUIJOTE-MFI data to assess the presence of AME in M31 following an independent procedure to that of B19 and Harper et al. (2023). The paper is organized as follows.

Section 2 describes the data used to build the M31 integrated spectral energy distribution (SED) and how they were processed. Section 3 describes the components assumed in the fitting of the SED of M31 and the fitting procedure itself. Section 4 presents the main results of this study, while Section 5 comments on the changes for these results with different assumptions. Finally, we present the main conclusions of this work in Section 6.

## 2 INPUT DATA

A summary of the maps used is presented in Table 1. All of them are used in HEALPix<fnref n="2" /> (Górski et al. 2005; Zonca et al. 2019) pixelization $N_{\text{side}} = 512$ and smoothed to a common resolution of 1°.

### 2.1 QUIJOTE-MFI data

The QUIJOTE CMB experiment (Rubiño-Martín et al. 2010) operates at the Teide Observatory (OT) of the Instituto de Astrofísica de Canarias (IAC), located at latitude $28^\circ 18'04''$ North and longitude $16^\circ 30'38''$ West. The first instrument installed on QUIJOTE was the Multi Frequency Instrument (MFI). This instrument observed the sky, in intensity and polarization, at 11, 13, 17, and 19 GHz, mostly in the so-called 'nominal mode' configuration, with continuous observations at constant elevation while spinning in azimuth. These observations constitute the Wide Survey (hereafter, WS), which is described in detail in Rubiño-Martín et al. (2023). In addition to the WS data set, we use new data obtained with raster mode observations specifically focused on M31, intended to improve the sensitivity in this area. All the analyses and results presented in this paper have been derived using 11 and 13 GHz maps. High noise levels in the 17 and 19 GHz intensity maps (see Table 2), mostly dictated by atmospheric 1/f noise, prevent the detection of emission from M31 at these frequencies so they are discarded for this analysis. The 11 and 13 GHz maps are produced combining the WS data with the data obtained in raster scan mode around the coordinates of M31. A single run of the QUIJOTE Map-Making code (Guidi et al. 2021) is required in order to optimize the recovering of the large angular-scale signal. The same strategy has been applied before (Guidi et al. 2023; Tramonte et al. 2023) and will be used in upcoming (González-González et al., in preparation; Ruiz-Granados et al., in preparation) publications using QUIJOTE-MFI data.

The data in raster-scan mode were taken between May and December 2016 for a total of 931 rasters accounting for 539.1 h of observations (see Table 3). These raster-scans are constant-elevation azimuth scans of width $\approx 12^\circ/\cos(\text{EL})$ and duration $\approx 35$ min, performed on different local coordinates while tracking the field, with elevations ranging between $32^\circ$ and $75^\circ$. After projection onto the sky plane, each of these observations results in a map of $\approx 12^\circ$ by $\approx 12^\circ$. The combination of all these observations results in a sky coverage as shown in the bottom-right-hand panel of Fig. 1, where the footprint of the raster scans is evident. This is one of the fields with highest integration time with the MFI. As a result, these are the maps with the best sensitivity of all of those obtained with the MFI to date, both in intensity and polarization. Final map sensitivities, computed through a jack-knife analysis on the Half Mission Difference Maps (HMDM; see Rubiño-Martín et al. 2023) are shown in Table 2: we can clearly see the improvement in sensitivity owing to the addition of raster data for most of the horns. In intensity (which are the data

<footnote n="2">²https://healpix.sourceforge.io</footnote>
<table><caption>**Table 1.** Summary of the surveys and frequency maps used in this analysis. The quoted photometry estimates were obtained as described in Section 3.2 after applying both the point sources and CMB subtractions described in Sections 3.4 and 3.5, respectively. The photometry estimates from B19, which are computed using the same apertures, are shown for comparison.</caption><thead><tr><th>Telescope</th><th>Frequency (GHz)</th><th>Calibration (per cent)</th><th>Resolution (arcmin)</th><th>S<sub>v</sub> from this work (Jy)</th><th>S<sub>v</sub> from Battistelli et al. (2019) (Jy)</th><th>Reference</th></tr></thead><tbody><tr><td>Various</td><td>0.408</td><td>10</td><td>51</td><td>13.7 ± 2.9</td><td>18.4 ± 1.6</td><td>Haslam et al. (1982)</td></tr><tr><td>Dwingeloo</td><td>0.82</td><td>10</td><td>72</td><td>8.6 ± 1.7</td><td>—</td><td>Remazeilles et al. (2015)</td></tr><tr><td>Stockert/Villa-Elisa</td><td>1.42</td><td>20</td><td>34.2</td><td>4.9 ± 1.3</td><td>5.28 ± 0.41</td><td>Berkhuijsen (1972)</td></tr><tr><td>QUIJOTE-MFI</td><td>11.2</td><td>5</td><td>53.2</td><td>1.79 ± 0.30</td><td>—</td><td>Reich (1982), Reich, Testori & Reich (2001)</td></tr><tr><td>QUIJOTE-MFI</td><td>12.9</td><td>5</td><td>53.5</td><td>1.94 ± 0.36</td><td>—</td><td>Rubiño-Martín et al. (2023)</td></tr><tr><td>WMAP</td><td>22.8</td><td>3</td><td>51.3</td><td>1.56 ± 0.14</td><td>2.00 ± 0.17</td><td>Bennett et al. (2013)</td></tr><tr><td>Planck</td><td>28.4</td><td>3</td><td>33.1</td><td>1.36 ± 0.14</td><td>1.86 ± 0.15</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>WMAP</td><td>33</td><td>3</td><td>39.1</td><td>1.26 ± 0.16</td><td>1.71 ± 0.21</td><td>Bennett et al. (2013)</td></tr><tr><td>WMAP</td><td>40.7</td><td>3</td><td>30.8</td><td>0.79 ± 0.15</td><td>1.31 ± 0.16</td><td>Bennett et al. (2013)</td></tr><tr><td>Planck</td><td>44.1</td><td>3</td><td>27.9</td><td>0.95 ± 0.18</td><td>1.45 ± 0.25</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>WMAP</td><td>60.7</td><td>3</td><td>21</td><td>0.84 ± 0.32</td><td>1.72 ± 0.42</td><td>Bennett et al. (2013)</td></tr><tr><td>Planck</td><td>70.4</td><td>3</td><td>13.1</td><td>1.01 ± 0.27</td><td>2.12 ± 0.36</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>WMAP</td><td>93.5</td><td>3</td><td>14.8</td><td>2.27 ± 0.82</td><td>3.5 ± 1.0</td><td>Bennett et al. (2013)</td></tr><tr><td>Planck</td><td>143</td><td>3</td><td>7.3</td><td>11.4 ± 1.1</td><td>15.7 ± 1.4</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>Planck</td><td>353</td><td>3</td><td>4.9</td><td>261 ± 22</td><td>318 ± 24</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>Planck</td><td>545</td><td>6.1</td><td>4.8</td><td>886 ± 83</td><td>1027 ± 73</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>Planck</td><td>857</td><td>6.4</td><td>4.6</td><td>2610 ± 240</td><td>3020 ± 190</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>COBE-DIRBE</td><td>1250</td><td>11.6</td><td>37.1</td><td>4680 ± 610</td><td>5330 ± 370 @ 1199 GHz</td><td>Hauser et al. (1998)</td></tr><tr><td>COBE-DIRBE</td><td>2143</td><td>10.6</td><td>38</td><td>6260 ± 730</td><td>7020 ± 230 @ 1874 GHz</td><td>Hauser et al. (1998)</td></tr><tr><td>COBE-DIRBE</td><td>3000</td><td>13.5</td><td>38.6</td><td>3010 ± 460</td><td>2980 ± 140 @ 2997 GHz</td><td>Hauser et al. (1998)</td></tr></tbody></table>

<table><caption>**Table 2.** Sensitivity estimates computed from a 3° radius circular region around M31 on HDMD, as is done in Rubiño-Martín et al. (2023). Estimates are in units of µK deg⁻¹. The estimates from the WS maps alone are also shown for comparison.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Data</th><th colspan="2">Horn 2</th><th colspan="2">Horn 3</th><th colspan="2">Horn 4</th></tr><tr><th>17 GHz</th><th>19 GHz</th><th>11 GHz</th><th>13 GHz</th><th>17 GHz</th><th>19 GHz</th></tr></thead><tbody><tr><td rowspan="2">I</td><td>WS+Raster</td><td>52.3</td><td>65.5</td><td>37.3</td><td>24.8</td><td>171.0</td><td>172.0</td></tr><tr><td>WS</td><td>173.4</td><td>204.8</td><td>96.9</td><td>69.0</td><td>217.8</td><td>237.2</td></tr><tr><td rowspan="2">Q</td><td>WS+Raster</td><td>41.2</td><td>54.1</td><td>13.4</td><td>11.1</td><td>12.4</td><td>10.8</td></tr><tr><td>WS</td><td>62.2</td><td>78.2</td><td>45.6</td><td>39.1</td><td>36.4</td><td>38.9</td></tr><tr><td rowspan="2">U</td><td>WS+Raster</td><td>42.0</td><td>52.9</td><td>13.9</td><td>11.3</td><td>11.8</td><td>11.1</td></tr><tr><td>WS</td><td>62.0</td><td>77.9</td><td>45.6</td><td>39.2</td><td>36.4</td><td>39.1</td></tr></tbody></table>

<table><caption>**Table 3.** Description of the QUIJOTE-MFI raster scans dedicated to M31. Third and fourth columns refer to the total amount of scans and time devoted to observe the M31 field, respectively. The fifth column shows the fraction of those scans (or time) that are merged together with the data from the WS. The last column shows the fraction of time coming from the WS (Rubiño-Martín et al. 2023) over the total after adding the raster data.</caption><thead><tr><th>Frequency (GHz)</th><th>Observing dates</th><th>Total number of scans</th><th>Total observing time (h)</th><th>Selected time fraction (per cent)</th><th>WS time fraction (per cent)</th></tr></thead><tbody><tr><td>11.2</td><td>May–December 2016</td><td>931</td><td>539.1</td><td>62.8</td><td>8.9</td></tr><tr><td>12.9</td><td></td><td></td><td></td><td>67.1</td><td>7.7</td></tr></tbody></table>

most of the results of this paper are based on), we reach sensitivities of 37 and 25 µK deg⁻¹, respectively at 11 and 13 GHz. These are better than those from the WS, 97 and 69 µK deg⁻¹, respectively. The improvement factor is then ≈2.7, slightly smaller than expected from the 10 times longer integration time. For the higher frequency horns, we can see that horn 2 behaves better in intensity while horn 4 does so in polarization, as was the case for the WS (Table 9, Rubiño-Martín et al. 2023).

We follow the same procedure and data processing pipeline that is currently being used in all QUIJOTE-MFI analyses dealing with raster observations (Guidi et al. 2023; Tramonte et al. 2023). Data are carefully inspected to identify and remove periods affected by systematic effects, in particular radio frequency interference (RFI)

and atmospheric contamination. With this aim, we produced Stokes IQU maps per horn and frequency of each individual observation and analysed their quality, first by eye and then by calculating their noise RMS. This allows us to identify particularly bad observations, which would degrade the final map quality. These observations are removed and not used to produce the stacked map. However, in the case of M31, this data cleaning on the maps was not sufficient, since those at 11 GHz in particular still presented clear signs of RFI, producing obvious striping along constant-declination directions. We then implemented an additional step consisting of combining the data in groups of azimuth and elevation. We then produced stacks of data for each group per horn and frequency by averaging the data in bins in azimuth. In these stacks, local signals clearly show
<table><caption>**Table 1.** Summary of the surveys and frequency maps used in this analysis. The quoted photometry estimates were obtained as described in Section 3.2 after applying both the point sources and CMB subtractions described in Sections 3.4 and 3.5, respectively. The photometry estimates from B19, which are computed using the same apertures, are shown for comparison.</caption><thead><tr><th>Telescope</th><th>Frequency (GHz)</th><th>Calibration (per cent)</th><th>Resolution (arcmin)</th><th>S<sub>v</sub> from this work (Jy)</th><th>S<sub>v</sub> from Battistelli et al. (2019) (Jy)</th><th>Reference</th></tr></thead><tbody><tr><td>Various</td><td>0.408</td><td>10</td><td>51</td><td>13.7 ± 2.9</td><td>18.4 ± 1.6</td><td>Haslam et al. (1982)</td></tr><tr><td>Dwingeloo</td><td>0.82</td><td>10</td><td>72</td><td>8.6 ± 1.7</td><td>—</td><td>Remazeilles et al. (2015)</td></tr><tr><td>Stockert/Villa-Elisa</td><td>1.42</td><td>20</td><td>34.2</td><td>4.9 ± 1.3</td><td>5.28 ± 0.41</td><td>Berkhuijsen (1972)</td></tr><tr><td>QUIJOTE-MFI</td><td>11.2</td><td>5</td><td>53.2</td><td>1.79 ± 0.30</td><td>—</td><td>Reich (1982), Reich, Testori & Reich (2001)</td></tr><tr><td>QUIJOTE-MFI</td><td>12.9</td><td>5</td><td>53.5</td><td>1.94 ± 0.36</td><td>—</td><td>Rubiño-Martín et al. (2023)</td></tr><tr><td>WMAP</td><td>22.8</td><td>3</td><td>51.3</td><td>1.56 ± 0.14</td><td>2.00 ± 0.17</td><td>Bennett et al. (2013)</td></tr><tr><td>Planck</td><td>28.4</td><td>3</td><td>33.1</td><td>1.36 ± 0.14</td><td>1.86 ± 0.15</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>WMAP</td><td>33</td><td>3</td><td>39.1</td><td>1.26 ± 0.16</td><td>1.71 ± 0.21</td><td>Bennett et al. (2013)</td></tr><tr><td>WMAP</td><td>40.7</td><td>3</td><td>30.8</td><td>0.79 ± 0.15</td><td>1.31 ± 0.16</td><td>Bennett et al. (2013)</td></tr><tr><td>Planck</td><td>44.1</td><td>3</td><td>27.9</td><td>0.95 ± 0.18</td><td>1.45 ± 0.25</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>WMAP</td><td>60.7</td><td>3</td><td>21</td><td>0.84 ± 0.32</td><td>1.72 ± 0.42</td><td>Bennett et al. (2013)</td></tr><tr><td>Planck</td><td>70.4</td><td>3</td><td>13.1</td><td>1.01 ± 0.27</td><td>2.12 ± 0.36</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>WMAP</td><td>93.5</td><td>3</td><td>14.8</td><td>2.27 ± 0.82</td><td>3.5 ± 1.0</td><td>Bennett et al. (2013)</td></tr><tr><td>Planck</td><td>143</td><td>3</td><td>7.3</td><td>11.4 ± 1.1</td><td>15.7 ± 1.4</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>Planck</td><td>353</td><td>3</td><td>4.9</td><td>261 ± 22</td><td>318 ± 24</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>Planck</td><td>545</td><td>6.1</td><td>4.8</td><td>886 ± 83</td><td>1027 ± 73</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>Planck</td><td>857</td><td>6.4</td><td>4.6</td><td>2610 ± 240</td><td>3020 ± 190</td><td>Planck Collaboration et al. (2020a)</td></tr><tr><td>COBE-DIRBE</td><td>1250</td><td>11.6</td><td>37.1</td><td>4680 ± 610</td><td>5330 ± 370 @ 1199 GHz</td><td>Hauser et al. (1998)</td></tr><tr><td>COBE-DIRBE</td><td>2143</td><td>10.6</td><td>38</td><td>6260 ± 730</td><td>7020 ± 230 @ 1874 GHz</td><td>Hauser et al. (1998)</td></tr><tr><td>COBE-DIRBE</td><td>3000</td><td>13.5</td><td>38.6</td><td>3010 ± 460</td><td>2980 ± 140 @ 2997 GHz</td><td>Hauser et al. (1998)</td></tr></tbody></table>

<table><caption>**Table 2.** Sensitivity estimates computed from a 3° radius circular region around M31 on HDMD, as is done in Rubiño-Martín et al. (2023). Estimates are in units of µK deg⁻¹. The estimates from the WS maps alone are also shown for comparison.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Data</th><th colspan="2">Horn 2</th><th colspan="2">Horn 3</th><th colspan="2">Horn 4</th></tr><tr><th>17 GHz</th><th>19 GHz</th><th>11 GHz</th><th>13 GHz</th><th>17 GHz</th><th>19 GHz</th></tr></thead><tbody><tr><td rowspan="2">I</td><td>WS+Raster</td><td>52.3</td><td>65.5</td><td>37.3</td><td>24.8</td><td>171.0</td><td>172.0</td></tr><tr><td>WS</td><td>173.4</td><td>204.8</td><td>96.9</td><td>69.0</td><td>217.8</td><td>237.2</td></tr><tr><td rowspan="2">Q</td><td>WS+Raster</td><td>41.2</td><td>54.1</td><td>13.4</td><td>11.1</td><td>12.4</td><td>10.8</td></tr><tr><td>WS</td><td>62.2</td><td>78.2</td><td>45.6</td><td>39.1</td><td>36.4</td><td>38.9</td></tr><tr><td rowspan="2">U</td><td>WS+Raster</td><td>42.0</td><td>52.9</td><td>13.9</td><td>11.3</td><td>11.8</td><td>11.1</td></tr><tr><td>WS</td><td>62.0</td><td>77.9</td><td>45.6</td><td>39.2</td><td>36.4</td><td>39.1</td></tr></tbody></table>

<table><caption>**Table 3.** Description of the QUIJOTE-MFI raster scans dedicated to M31. Third and fourth columns refer to the total amount of scans and time devoted to observe the M31 field, respectively. The fifth column shows the fraction of those scans (or time) that are merged together with the data from the WS. The last column shows the fraction of time coming from the WS (Rubiño-Martín et al. 2023) over the total after adding the raster data.</caption><thead><tr><th>Frequency (GHz)</th><th>Observing dates</th><th>Total number of scans</th><th>Total observing time (h)</th><th>Selected time fraction (per cent)</th><th>WS time fraction (per cent)</th></tr></thead><tbody><tr><td>11.2</td><td>May–December 2016</td><td>931</td><td>539.1</td><td>62.8</td><td>8.9</td></tr><tr><td>12.9</td><td></td><td></td><td></td><td>67.1</td><td>7.7</td></tr></tbody></table>

most of the results of this paper are based on), we reach sensitivities of 37 and 25 µK deg⁻¹, respectively at 11 and 13 GHz. These are better than those from the WS, 97 and 69 µK deg⁻¹, respectively. The improvement factor is then ≈2.7, slightly smaller than expected from the 10 times longer integration time. For the higher frequency horns, we can see that horn 2 behaves better in intensity while horn 4 does so in polarization, as was the case for the WS (Table 9, Rubiño-Martín et al. 2023).

We follow the same procedure and data processing pipeline that is currently being used in all QUIJOTE-MFI analyses dealing with raster observations (Guidi et al. 2023; Tramonte et al. 2023). Data are carefully inspected to identify and remove periods affected by systematic effects, in particular radio frequency interference (RFI)

and atmospheric contamination. With this aim, we produced Stokes IQU maps per horn and frequency of each individual observation and analysed their quality, first by eye and then by calculating their noise RMS. This allows us to identify particularly bad observations, which would degrade the final map quality. These observations are removed and not used to produce the stacked map. However, in the case of M31, this data cleaning on the maps was not sufficient, since those at 11 GHz in particular still presented clear signs of RFI, producing obvious striping along constant-declination directions. We then implemented an additional step consisting of combining the data in groups of azimuth and elevation. We then produced stacks of data for each group per horn and frequency by averaging the data in bins in azimuth. In these stacks, local signals clearly show
<figure><img src="image_1.png" /><figcaption>**Figure 1.** QUIJOTE-MFI 11 GHz intensity (top row) and $N_{\text{hits}}$ (bottom row) maps for the WS release (left column), and WS plus raster scan data on the M31 field (right column). The hits are defined as 40 ms samples in a HEALPix map of $N_{\text{side}} = 512$. The new raster data account for an exposure time approximately 10 times longer than that from the WS. The elliptical and rectangular regions in black are those used to compute M31 density fluxes (and their uncertainties) shown in Table 1 and used in Section 4. The photometric results in these maps using those regions are also quoted in the top row maps. In the bottom row maps, we can see how the exposure time increases by an order of magnitude when including the raster data with that from the WS.</figcaption></figure>

up as bumps at specific AZ ranges. A template could then be built and removed from each scan, but we preferred to follow a more conservative approach and instead removed specific azimuth ranges of horns and frequencies that presented signs of RFI. On top of this, all observations performed at AZ<180° and EL ≈ 32° were removed as they presented particularly strong RFI that was very difficult to correct. All these processes resulted in removing 37 per cent of the 11 GHz data and 33 per cent of the 13 GHz data, as indicated in Table 3. Even after this, raster data account for more than 90 per cent of the data considered in the final maps.

## 2.2 WMAP/Planck/DIRBE

We used the Wilkinson Microwave Anisotropy Probe (WMAP; Bennett et al. 2013) 9-year and Planck PR3 2018 maps (Planck Collaboration et al. 2020a). We also used the zodiacal light subtracted version of COBE Diffuse Infrared Background Experiment (COBE-DIRBE; Hauser et al. 1998) maps at 240, 140, and 100 µm. WMAP and COBE-DIRBE data are available in the Legacy Archive for Microwave Background Data Analysis (LAMBDA<fnref n="3" />), while Planck data are available through the Planck Legacy Archive (PLA<fnref n="4" />). The calibration uncertainties for WMAP and those Planck bands calibrated against the CMB dipole (up to 353 GHz included) are

increased to 3 per cent.<fnref n="5" /> Note that, nominally, WMAP and Planck data have much better (sub-per cent) uncertainties. However, we decided to use more conservative values to account not only for the global gain calibration uncertainties but also for issues related with beam inaccuracies or colour corrections. These values are the same as those used in recent QUIJOTE papers (e.g. López-Caraballo et al., 2023; Poidevin et al. 2023; Tramonte et al. 2023) or in other papers (e.g. Planck Collaboration et al. 2014a). Planck maps at 100 and 217 GHz are not used because of contamination by CO residuals.

## 2.3 Other ancillary surveys

We also used three low-frequency surveys covering the M31 region: the Haslam et al. (1982) map at 0.408 GHz, the Berkhuijsen (1972) map at 0.820 GHz and the Reich (1982), Reich & Reich (1988) map at 1.42 GHz. We increased the calibration uncertainty from Berkhuijsen (1972) from 6 to 10 per cent to account for the fact that its parent angular resolution (72 arcmin) is slightly larger than the 1° resolution used in this analysis. Besides, the Reich (1982) map requires applying a recalibration factor when studying structures on the main beam scale. This is normally assumed to be 1.55 (e.g. Reich & Reich 1988) when studying point sources, but decreases when focusing on resolved objects such as M31 (Fernández-Torreiro

<footnote n="3">³https://lambda.gsfc.nasa.gov</footnote>

<footnote n="4">⁴http://pla.esac.esa.int/pla/</footnote>

<footnote n="5">⁵Planck -HFI 545 and 857 GHz bands, which are calibrated using planets, have calibration uncertainties greater than 6 per cent.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** QUIJOTE-MFI 11 GHz intensity (top row) and $N_{\text{hits}}$ (bottom row) maps for the WS release (left column), and WS plus raster scan data on the M31 field (right column). The hits are defined as 40 ms samples in a HEALPix map of $N_{\text{side}} = 512$. The new raster data account for an exposure time approximately 10 times longer than that from the WS. The elliptical and rectangular regions in black are those used to compute M31 density fluxes (and their uncertainties) shown in Table 1 and used in Section 4. The photometric results in these maps using those regions are also quoted in the top row maps. In the bottom row maps, we can see how the exposure time increases by an order of magnitude when including the raster data with that from the WS.</figcaption></figure>

up as bumps at specific AZ ranges. A template could then be built and removed from each scan, but we preferred to follow a more conservative approach and instead removed specific azimuth ranges of horns and frequencies that presented signs of RFI. On top of this, all observations performed at AZ<180° and EL ≈ 32° were removed as they presented particularly strong RFI that was very difficult to correct. All these processes resulted in removing 37 per cent of the 11 GHz data and 33 per cent of the 13 GHz data, as indicated in Table 3. Even after this, raster data account for more than 90 per cent of the data considered in the final maps.

## 2.2 WMAP/Planck/DIRBE

We used the Wilkinson Microwave Anisotropy Probe (WMAP; Bennett et al. 2013) 9-year and Planck PR3 2018 maps (Planck Collaboration et al. 2020a). We also used the zodiacal light subtracted version of COBE Diffuse Infrared Background Experiment (COBE-DIRBE; Hauser et al. 1998) maps at 240, 140, and 100 µm. WMAP and COBE-DIRBE data are available in the Legacy Archive for Microwave Background Data Analysis (LAMBDA<fnref n="3" />), while Planck data are available through the Planck Legacy Archive (PLA<fnref n="4" />). The calibration uncertainties for WMAP and those Planck bands calibrated against the CMB dipole (up to 353 GHz included) are

increased to 3 per cent.<fnref n="5" /> Note that, nominally, WMAP and Planck data have much better (sub-per cent) uncertainties. However, we decided to use more conservative values to account not only for the global gain calibration uncertainties but also for issues related with beam inaccuracies or colour corrections. These values are the same as those used in recent QUIJOTE papers (e.g. López-Caraballo et al., 2023; Poidevin et al. 2023; Tramonte et al. 2023) or in other papers (e.g. Planck Collaboration et al. 2014a). Planck maps at 100 and 217 GHz are not used because of contamination by CO residuals.

## 2.3 Other ancillary surveys

We also used three low-frequency surveys covering the M31 region: the Haslam et al. (1982) map at 0.408 GHz, the Berkhuijsen (1972) map at 0.820 GHz and the Reich (1982), Reich & Reich (1988) map at 1.42 GHz. We increased the calibration uncertainty from Berkhuijsen (1972) from 6 to 10 per cent to account for the fact that its parent angular resolution (72 arcmin) is slightly larger than the 1° resolution used in this analysis. Besides, the Reich (1982) map requires applying a recalibration factor when studying structures on the main beam scale. This is normally assumed to be 1.55 (e.g. Reich & Reich 1988) when studying point sources, but decreases when focusing on resolved objects such as M31 (Fernández-Torreiro

<footnote n="3">³https://lambda.gsfc.nasa.gov</footnote>

<footnote n="4">⁴http://pla.esac.esa.int/pla/</footnote>

<footnote n="5">⁵Planck -HFI 545 and 857 GHz bands, which are calibrated using planets, have calibration uncertainties greater than 6 per cent.</footnote>
et al. 2023). Therefore, we kept the 1.3 value used by Planck Collaboration et al. (2015b) and increased its calibration uncertainty to 20 per cent. We used the Remazeilles et al. (2015) and Paradis et al. (2012) reprocessed versions of the Haslam et al. (1982) and Reich (1982); Reich & Reich (1988) maps, respectively. The three maps described in this section are available in the LAMBDA repository.

## 2.4 High resolution surveys

Several telescopes have observed M31 with arcminute resolution in the past. A review of these observations is presented in Berkhuijsen, Beck & Hoernes (2003). However, the fields covered are too small to be introduced into the analysis. The same happens for SRT data from B19 and Fatigoni et al. (2021): when smoothing the data to 1° resolutions, artefacts due to the presence of edges in the image may arise, which further complicate the analysis. Because of this, we are not using any of these higher resolution data in our analysis and the final set of maps used to constrain the SEDs described in Section 3 is the one from Table 1. We provide only qualitative comparisons to these higher-resolution analyses.

# 3 METHODOLOGY

## 3.1 SED modelling

We have considered up to four different emission components in this study: synchrotron, free-free, AME and thermal dust. These are modelled as in Fernández-Torreiro et al. (2023), with the main difference being the lack of a component accounting for CMB anisotropies. Instead, we subtract the CMB anisotropies using the NILC (see Section 3.5, Planck Collaboration et al. 2020b) map. Therefore, the flux density at frequency $v$ is described as:

$$ 
\begin{aligned}
S_{\nu}^{\text{total}}(\theta) &= S_{\nu}^{\text{syn}}(A_{\text{1 GHz}}, \alpha_{\text{syn}}) + S_{\nu}^{\text{ff}}(\text{EM}) \\
&\quad + S_{\nu}^{\text{AME}}(A_{\text{AME}}, v_{\text{AME}}, W_{\text{AME}}) \\
&\quad + S_{\nu}^{\text{dust}}(\tau_{\text{353}}, \beta_{\text{d}}, T_{\text{d}}) \\
&= A_{\text{1 GHz}} \left(\frac{\nu}{1 \text{ GHz}}\right)^{\alpha_{\text{syn}}} + \frac{2k_{\text{B}}\nu^2}{c^2} \Omega T_{\text{ff}} \\
&\quad + A_{\text{AME}} \exp \left[ -\frac{1}{2W_{\text{AME}}^2} \ln^2 \left(\frac{\nu}{v_{\text{AME}}}\right) \right] \\
&\quad + \frac{2\nu h v^3}{c^2} \left(\frac{\nu}{353 \text{ GHz}}\right)^{\beta_d} \tau_{\text{353}} (e^{\nu h/k_{\text{B}} T_d} - 1)^{-1} \Omega,
\end{aligned}
\quad (1) $$

where $A_{1\text{GHz}}$ stands for the synchrotron flux density at 1 GHz and $\alpha_{\text{syn}}$ for the synchrotron spectral index; EM for the emission measure, which is used in the definition of $T_{\text{ff}}$; $A_{\text{AME}}$, $v_{\text{AME}}$, and $W_{\text{AME}}$ for the AME peak flux density, frequency, and width; and $\tau_{\text{353}}$, $\beta_{\text{d}}$, and $T_{\text{d}}$ for the opacity, the spectral index, and the dust grain temperature, respectively, which describe the spectrum of the thermal dust emission. $\Omega$ accounts for the solid angle of the M31 aperture (which is the same used in B19): $\Omega = 4.38 \text{ deg}^2$.

Note that, as it has now become the standard in similar analyses, for the AME we resort to the phenomenological model proposed by Stevenson (2014). This model provides a good fit to most of the spinning-dust models. For the thermal dust emission, we use a single-temperature modified black body. We are aware that in the integrated SED of M31, we mix different dust components, but even in this case this model provides a reasonably good fit. In any case, a multicomponent fit would not be possible given the sparse spectral coverage.

## 3.2 Aperture photometry

We adopt the same M31 aperture used in B19: an ellipse with semi-axes 91.5 and 59.5 arcmin, position angle $-52^\circ$ (east-to-north) and centred on (RAJ2000, DecJ2000) = (10.68°, 41.27°).<fnref n="6" /> Our background aperture is defined as those pixels also present in the B19 SRT image outside the M31 ellipse; this defines a rectangle with edges located at RAJ2000 = (9.15°, 12.25°) and DECJ2000 = (39.74°, 42.82°). We decided to keep the source and background regions identical to those in B19 not only to provide a direct comparison with that study, but also to be able to use the point source catalogue released later by Fatigoni et al. (2021). Both M31 and background apertures can be seen in Fig. 2. We computed M31 flux densities as in previous studies (e.g. Rubiño-Martín et al. 2012; Génova-Santos et al. 2017; Fernández-Torreiro et al. 2023) by using the mean and the median values within the source and background apertures, respectively [equation (2)]. The uncertainty for these flux densities is computed from the standard deviation within the background aperture [equation (3)]:

$$ S_{\nu} = a(\nu)\Omega T = \frac{2k_{\text{B}}\nu^2}{c^2} \frac{x^2 e^x}{(e^x - 1)^2} \Omega [\bar{T}_{\text{aper}} - \text{med}(T_{\text{BG}})], \quad (2) $$

$$ \sigma_{\text{AP}} = a(\nu)\Omega \sigma(T_{\text{BG}}) \sqrt{\frac{n_{\text{beam}}}{n_{\text{aper}}} + \frac{\pi n_{\text{beam}}}{2 n_{\text{BG}}}}, \quad (3) $$

where $n_{\text{beam}}$ is the number of pixels within a 1° beam, while $n_{\text{aper}}$ and $n_{\text{BG}}$ are the number of pixels within the M31 and background apertures, respectively. We add this factor because of the assumption that the noise is completely correlated on beam scales: this is a conservative way of estimating the uncertainties for aperture photometry studies (Génova-Santos et al. 2015, 2017). Finally, the calibration uncertainty quoted in Table 1 is added quadratically to the uncertainty estimation:

$$ \sigma_{S_{\nu}} = \sqrt{\sigma_{\text{AP}}^2 + \text{cal}^2 \cdot S_{\nu}^2}. \quad (4) $$

## 3.3 MCMC

We ran several maximum likelihood estimator (MLE) analyses using Markov Chain Monte Carlo (MCMC) samplers for various scenarios, using emcee (Foreman-Mackey et al. 2013). Each of these MCMCs ran for a fixed length of $N = 10^5$ steps. Once they ended, the convergence of the chains was assessed by computing their autocorrelation times, $\tau$. The chains have converged if $\tau < N/50$: all cases fulfilled this requirement. Colour corrections (cc) are embedded within the MCMC log-likelihood ($\log L$) and computed iteratively, using FASTCC and INTERPCC from Peel et al. (2022), as in the following equation:

$$ \log L = -0.5 \cdot \left( S - \frac{S^{\text{total}}(\theta)}{cc} \right)^{\text{T}} C^{-1} \left( S - \frac{S^{\text{total}}(\theta)}{cc} \right), \quad (5) $$

where $S$ and $S^{\text{total}}(\theta)$ stand for the measured and expected flux densities. $C$ is the covariance matrix between the surveys, which has all its off-diagonal elements equal to zero except for those corresponding to the 11 and 13 GHz maps, whose noise is partially correlated. Further details on this methodology are available in Section 3.3.2 of Fernández-Torreiro et al. (2023). We use the same flat minimal priors on the parameter posteriors as in that study.

<footnote n="6">⁶https://simbad.u-strasbg.fr/simbad/sim-basic?Ident=M31&submit = SIMBAD + search</footnote>
et al. 2023). Therefore, we kept the 1.3 value used by Planck Collaboration et al. (2015b) and increased its calibration uncertainty to 20 per cent. We used the Remazeilles et al. (2015) and Paradis et al. (2012) reprocessed versions of the Haslam et al. (1982) and Reich (1982); Reich & Reich (1988) maps, respectively. The three maps described in this section are available in the LAMBDA repository.

## 2.4 High resolution surveys

Several telescopes have observed M31 with arcminute resolution in the past. A review of these observations is presented in Berkhuijsen, Beck & Hoernes (2003). However, the fields covered are too small to be introduced into the analysis. The same happens for SRT data from B19 and Fatigoni et al. (2021): when smoothing the data to 1° resolutions, artefacts due to the presence of edges in the image may arise, which further complicate the analysis. Because of this, we are not using any of these higher resolution data in our analysis and the final set of maps used to constrain the SEDs described in Section 3 is the one from Table 1. We provide only qualitative comparisons to these higher-resolution analyses.

# 3 METHODOLOGY

## 3.1 SED modelling

We have considered up to four different emission components in this study: synchrotron, free-free, AME and thermal dust. These are modelled as in Fernández-Torreiro et al. (2023), with the main difference being the lack of a component accounting for CMB anisotropies. Instead, we subtract the CMB anisotropies using the NILC (see Section 3.5, Planck Collaboration et al. 2020b) map. Therefore, the flux density at frequency $v$ is described as:

$$ 
\begin{aligned}
S_{\nu}^{\text{total}}(\theta) &= S_{\nu}^{\text{syn}}(A_{1\,\text{GHz}}, \alpha_{\text{syn}}) + S_{\nu}^{\text{ff}}(\text{EM}) \\
&\quad + S_{\nu}^{\text{AME}}(A_{\text{AME}}, v_{\text{AME}}, W_{\text{AME}}) \\
&\quad + S_{\nu}^{\text{dust}}(\tau_{353}, \beta_{d}, T_{d}) \\
&= A_{1\,\text{GHz}} \left(\frac{\nu}{1\,\text{GHz}}\right)^{\alpha_{\text{syn}}} + \frac{2k_{\text{B}}\nu^2}{c^2} \Omega T_{\text{ff}} \\
&\quad + A_{\text{AME}} \exp \left[ -\frac{1}{2W_{\text{AME}}^2} \ln^2 \left(\frac{\nu}{v_{\text{AME}}}\right) \right] \\
&\quad + \frac{2\nu h v^3}{c^2} \left(\frac{\nu}{353\,\text{GHz}}\right)^{\beta_d} \tau_{353} (e^{hv/k_{\text{B}}T_d} - 1)^{-1} \Omega,
\end{aligned}
\quad (1) $$

where $A_{1\,\text{GHz}}$ stands for the synchrotron flux density at 1 GHz and $\alpha_{\text{syn}}$ for the synchrotron spectral index; EM for the emission measure, which is used in the definition of $T_{\text{ff}}$; $A_{\text{AME}}$, $v_{\text{AME}}$, and $W_{\text{AME}}$ for the AME peak flux density, frequency, and width; and $\tau_{353}$, $\beta_d$, and $T_d$ for the opacity, the spectral index, and the dust grain temperature, respectively, which describe the spectrum of the thermal dust emission. $\Omega$ accounts for the solid angle of the M31 aperture (which is the same used in B19): $\Omega = 4.38 \, \text{deg}^2$.

Note that, as it has now become the standard in similar analyses, for the AME we resort to the phenomenological model proposed by Stevenson (2014). This model provides a good fit to most of the spinning-dust models. For the thermal dust emission, we use a single-temperature modified black body. We are aware that in the integrated SED of M31, we mix different dust components, but even in this case this model provides a reasonably good fit. In any case, a multicomponent fit would not be possible given the sparse spectral coverage.

## 3.2 Aperture photometry

We adopt the same M31 aperture used in B19: an ellipse with semi-axes 91.5 and 59.5 arcmin, position angle $-52^\circ$ (east-to-north) and centred on (RAJ2000, DecJ2000) = (10.68$^\circ$, 41.27$^\circ)$.<fnref n="6" /> Our background aperture is defined as those pixels also present in the B19 SRT image outside the M31 ellipse; this defines a rectangle with edges located at RAJ2000 = (9.15$^\circ$, 12.25$^\circ)$ and DECJ2000 = (39.74$^\circ$, 42.82$^\circ)$. We decided to keep the source and background regions identical to those in B19 not only to provide a direct comparison with that study, but also to be able to use the point source catalogue released later by Fatigoni et al. (2021). Both M31 and background apertures can be seen in Fig. 2. We computed M31 flux densities as in previous studies (e.g. Rubiño-Martín et al. 2012; Génova-Santos et al. 2017; Fernández-Torreiro et al. 2023) by using the mean and the median values within the source and background apertures, respectively [equation (2)]. The uncertainty for these flux densities is computed from the standard deviation within the background aperture [equation (3)]:

$$ S_{\nu} = a(\nu)\Omega T = \frac{2k_{\text{B}}\nu^2}{c^2} \frac{x^2 e^x}{(e^x - 1)^2} \Omega [\bar{T}_{\text{aper}} - \text{med}(T_{\text{BG}})], \quad (2) $$

$$ \sigma_{\text{AP}} = a(\nu)\Omega \sigma(T_{\text{BG}}) \sqrt{\frac{n_{\text{beam}}}{n_{\text{aper}}} + \frac{\pi n_{\text{beam}}}{2 n_{\text{BG}}}}, \quad (3) $$

where $n_{\text{beam}}$ is the number of pixels within a 1° beam, while $n_{\text{aper}}$ and $n_{\text{BG}}$ are the number of pixels within the M31 and background apertures, respectively. We add this factor because of the assumption that the noise is completely correlated on beam scales: this is a conservative way of estimating the uncertainties for aperture photometry studies (Génova-Santos et al. 2015, 2017). Finally, the calibration uncertainty quoted in Table 1 is added quadratically to the uncertainty estimation:

$$ \sigma_{S_{\nu}} = \sqrt{\sigma_{\text{AP}}^2 + \text{cal}^2 \cdot S_{\nu}^2}. \quad (4) $$

## 3.3 MCMC

We ran several maximum likelihood estimator (MLE) analyses using Markov Chain Monte Carlo (MCMC) samplers for various scenarios, using emcee (Foreman-Mackey et al. 2013). Each of these MCMCs ran for a fixed length of $N = 10^5$ steps. Once they ended, the convergence of the chains was assessed by computing their autocorrelation times, $\tau$. The chains have converged if $\tau < N/50$: all cases fulfilled this requirement. Colour corrections (cc) are embedded within the MCMC log-likelihood ($\log L$) and computed iteratively, using FASTCC and INTERPCC from Peel et al. (2022), as in the following equation:

$$ \log L = -0.5 \cdot \left( \mathbf{S} - \frac{\mathbf{S}^{\text{total}}(\theta)}{\mathbf{cc}} \right)^{\mathrm{T}} \mathbf{C}^{-1} \left( \mathbf{S} - \frac{\mathbf{S}^{\text{total}}(\theta)}{\mathbf{cc}} \right), \quad (5) $$

where $\mathbf{S}$ and $\mathbf{S}^{\text{total}}(\theta)$ stand for the measured and expected flux densities. $\mathbf{C}$ is the covariance matrix between the surveys, which has all its off-diagonal elements equal to zero except for those corresponding to the 11 and 13 GHz maps, whose noise is partially correlated. Further details on this methodology are available in Section 3.3.2 of Fernández-Torreiro et al. (2023). We use the same flat minimal priors on the parameter posteriors as in that study.

<footnote n="6">⁶https://simbad.u-strasbg.fr/simbad/sim-basic?Ident=M31&submit = SIMBAD + search</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Examples of maps before and after CMB and point source subtraction. The first column shows the original maps, while the second shows the maps after subtracting the CMB. The third column shows the emission from point sources at the frequencies of the maps. The last column shows the final processed maps, which are used to compute the flux densities. The different rows host different maps: from top to bottom, those are Reich at 1.42 GHz, QUIJOTE-MFI at 11.2 and 12.9 GHz, WMAP at 22.8 GHz, and Planck-HFI at 143 GHz. The elliptical aperture used to extract M31 flux is shown as a solid line. The rectangular region marked by the dashed line is used to compute the background emission level. All maps have mK units. The colourbar levels for the point source maps are fixed at 25 per cent of the amplitude of the colourbar in the original maps. In WMAP K and MFI 11 GHz maps, it is clear how 5C+3.50 (Section 3.4.1) emission distorts the elliptical morphology of M31 as seen in SRT (Fig. 3) or thermal dust (143 GHz, bottom row) maps. This morphology is recovered after subtracting the source, as can be seen in the last column. The most important correction for the 143 GHz map, however, is the CMB anisotropy subtraction.</figcaption></figure>

## 3.4 Point source subtraction

We used the catalogue from Fatigoni et al. (2021)<fnref n="7" /> to subtract the point sources within the SRT M31 observed field (Fig. 3). The catalogue provides the flux density at $v_{ref} = 1$ GHz (A) plus a constant spectral index ($\alpha$) and a curvature parameter ($k$), so the flux density at a frequency $v$ is defined as:

$$S = A \left( \frac{v}{v_{ref}} \right)^{\alpha} \exp \left[ k \left( \frac{v}{v_{ref}} \right)^{-1} \right]. \qquad (6)$$

### 3.4.1 B3 0035+413 variability

One of the point sources, B3 0035+413 (or 5C+3.50), required a detailed analysis for two reasons. First, its spectral index is one of the flattest in the full catalogue, with $\alpha = -0.08$,<fnref n="8" /> compared

to the median value $\alpha = -0.78$ across all sources. This, together with its high flux density, implies that this source remains bright for many of the frequency bands studied. It has a null curvature ($k = 0$) estimate. It is clear in Fig. 2 how the morphology of M31 in the WMAP K map resembles the ellipse seen in SRT or thermal dust data only once the emission from this source is subtracted. Second, this is a variable source: excess emission was recorded by the 40-m telescope at the Owens Valley Radio Observatory (OVRO) at 15 GHz between 2009 and 2012 (Richards et al. 2011), increasing from 0.45 to 0.55 Jy. This excess overlaps in time with Planck observations and therefore accounts for 20 per cent of its flux before the event. It is impossible to know the exact behaviour of this source over the

frequencies. Variable extragalactic sources tend to have flatter spectra (Rys & Machalski 1990; Healey et al. 2007) because of optically thick synchrotron emission. This requires high surface brightness, which for the range of flux densities that we actually measure, implies small angular sizes. On the other hand, variability on yearly time-scales require small light-crossing time-scales, or sizes. A source with a steep spectrum has more flux coming from its outer regions, thus its light-crossing time-scale would be large and variability could not be measured.

<footnote n="7">⁷https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/651/A98</footnote>

<footnote n="8">⁸No uncertainties are provided in the catalogue. Harper et al. (2023) also found a significantly flat spectrum ($\alpha = -0.04 \pm 0.01$) at radio-to-microwave</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Examples of maps before and after CMB and point source subtraction. The first column shows the original maps, while the second shows the maps after subtracting the CMB. The third column shows the emission from point sources at the frequencies of the maps. The last column shows the final processed maps, which are used to compute the flux densities. The different rows host different maps: from top to bottom, those are Reich at 1.42 GHz, QUIJOTE-MFI at 11.2 and 12.9 GHz, WMAP at 22.8 GHz, and Planck-HFI at 143 GHz. The elliptical aperture used to extract M31 flux is shown as a solid line. The rectangular region marked by the dashed line is used to compute the background emission level. All maps have mK units. The colourbar levels for the point source maps are fixed at 25 per cent of the amplitude of the colourbar in the original maps. In WMAP K and MFI 11 GHz maps, it is clear how 5C+3.50 (Section 3.4.1) emission distorts the elliptical morphology of M31 as seen in SRT (Fig. 3) or thermal dust (143 GHz, bottom row) maps. This morphology is recovered after subtracting the source, as can be seen in the last column. The most important correction for the 143 GHz map, however, is the CMB anisotropy subtraction.</figcaption></figure>

## 3.4 Point source subtraction

We used the catalogue from Fatigoni et al. (2021)<fnref n="7" /> to subtract the point sources within the SRT M31 observed field (Fig. 3). The catalogue provides the flux density at $v_{ref} = 1$ GHz (A) plus a constant spectral index ($\alpha$) and a curvature parameter ($k$), so the flux density at a frequency $v$ is defined as:

$$S = A \left( \frac{v}{v_{ref}} \right)^{\alpha} \exp \left[ k \left( \frac{v}{v_{ref}} \right)^{-1} \right]. \qquad (6)$$

### 3.4.1 B3 0035+413 variability

One of the point sources, B3 0035+413 (or 5C+3.50), required a detailed analysis for two reasons. First, its spectral index is one of the flattest in the full catalogue, with $\alpha = -0.08$,<fnref n="8" /> compared

to the median value $\alpha = -0.78$ across all sources. This, together with its high flux density, implies that this source remains bright for many of the frequency bands studied. It has a null curvature ($k = 0$) estimate. It is clear in Fig. 2 how the morphology of M31 in the WMAP K map resembles the ellipse seen in SRT or thermal dust data only once the emission from this source is subtracted. Second, this is a variable source: excess emission was recorded by the 40-m telescope at the Owens Valley Radio Observatory (OVRO) at 15 GHz between 2009 and 2012 (Richards et al. 2011), increasing from 0.45 to 0.55 Jy. This excess overlaps in time with Planck observations and therefore accounts for 20 per cent of its flux before the event. It is impossible to know the exact behaviour of this source over the

frequencies. Variable extragalactic sources tend to have flatter spectra (Rys & Machalski 1990; Healey et al. 2007) because of optically thick synchrotron emission. This requires high surface brightness, which for the range of flux densities that we actually measure, implies small angular sizes. On the other hand, variability on yearly time-scales require small light-crossing time-scales, or sizes. A source with a steep spectrum has more flux coming from its outer regions, thus its light-crossing time-scale would be large and variability could not be measured.

<footnote n="7">⁷https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/651/A98</footnote>

<footnote n="8">⁸No uncertainties are provided in the catalogue. Harper et al. (2023) also found a significantly flat spectrum ($\alpha = -0.04 \pm 0.01$) at radio-to-microwave</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Left: M31 map obtained with the SRT at 6.6 GHz before subtracting point sources. Point sources are indicated with black circles. Middle: point sources map at 6.6 GHz SRT resolution (2.9 arcmin), according to Fatigoni et al. (2021). Right: M31 map after subtracting the point sources. In this last map, the opacity for the point source markers is directly proportional to their amplitude at 6.6 GHz.</figcaption></figure>

whole time span covered by the surveys depicted in Table 1 (from the '80s well into the 2010s). Besides, variations can be limited to only part of the source spectra. Therefore, we decided to estimate the change in the photometric results induced by large variations (worst case scenario) of this source flux, and increase the uncertainties accordingly.

We ran a set of $10^4$ simulations for this source, defining its amplitude at 6.6 GHz from a random Gaussian realization assuming a $\sigma$ equal to 50 per cent of the measured flux in the Fatigoni et al. (2021) catalogue. We assumed this variability level after the analysis ran in Section 4.2 of Harper et al. (2023), where the variability of this source was found to be between 30 and 50 per cent from yearly Planck 28.4 GHz maps. We used the 5C+3.50 broken power spectral model from Harper et al. (2023) to build an estimate for this source intensity at each map frequency. This estimate is then projected onto a map for each frequency assuming the SRT FWHM beam. These maps are smoothed to 1° resolution and transformed to HEALPix maps. These maps, which contain only the signal from 5C+3.50, are then subtracted from the frequency maps. We performed aperture photometry on these final maps, as explained in Section 3.2, and assessed the changes in the extracted flux density for the full M31 aperture. We found that the introduced 50 per cent $\sigma$ variability level translates into a maximum 40 per cent variation in the computed flux density of M31 in the Planck 28.4 GHz band. We increased the flux density error at each frequency by adding in quadrature the uncertainty associated with the variability of this source. Therefore, equation (4) transforms into:

$$ \sigma_{S_v} = \sqrt{\sigma_{\text{AP}}^2 + \text{cal}^2 \cdot S_v^2 + (\sigma_v^{\text{var}})^2 \cdot S_v^2}, \quad (7) $$

where $\sigma_v^{\text{var}}$ is the variation level found for each band. It is only larger than 10 per cent for WMAP and Planck bands between 20 and 50 GHz, while being negligible for low ($\nu < 10$ GHz) and high ($\nu > 90$ GHz) frequency bands.

## 3.5 Subtraction of CMB anisotropies

The contribution of CMB anisotropies to the global SED of M31 is not negligible, as can be seen in the Planck 143 GHz maps in the bottom row of Fig. 2. In fact, given the relatively large angular size of M31 and its low flux density, CMB anisotropies make an important contribution to the uncertainty in the measured flux densities. As mentioned in Section 2, we used 1° smoothed maps in this study: at that scale, CMB anisotropies show an RMS of $\approx 80$ µK. When translated into flux density, that uncertainty value dominates over the dispersion of the rest of the emission from the background region. It is thus preferable first to subtract the CMB anisotropies and then fit for the other components at the SED level. Several maps tracing CMB anisotropies exist, but we focused on the four official Planck collaboration CMB maps from SMICA, NILC, SEVEM, and COMMANDER methods (Planck Collaboration et al. 2020b).

In order to study to what degree the results presented in this paper could be affected by the selection of the CMB component-separated map, we performed aperture photometry on all CMB maps. First, we used SMICA, following the methodology of B19. We found that the SMICA CMB integrated contribution in the M31 region was the most important of the four maps (+1.0 µK). This estimate (and those for the other maps quoted later on) was computed as the mean temperature within the M31 aperture minus the median of the background region, using the apertures from Section 3.2. COMMANDER is the opposite case: its estimate is the least important of the four (-2.4 µK). NILC and SEVEM yielded similar CMB mean temperatures (-1.3 µK and -0.8 µK), which sit between the SMICA and COMMANDER values, as good trade-off values. We decided to use NILC after performing a visual inspection of the maps (there are visible residuals in SEVEM on the M31 outer ring). Using any of the extreme values (SMICA or COMMANDER) biases the AME estimation to high or low values. This is because the amplitude from the CMB maps within the M31 aperture is comparable to the one from the emission coming from the galaxy itself. This is further discussed in Section 5.1. Finally, in order to
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: M31 map obtained with the SRT at 6.6 GHz before subtracting point sources. Point sources are indicated with black circles. Middle: point sources map at 6.6 GHz SRT resolution (2.9 arcmin), according to Fatigoni et al. (2021). Right: M31 map after subtracting the point sources. In this last map, the opacity for the point source markers is directly proportional to their amplitude at 6.6 GHz.</figcaption></figure>

whole time span covered by the surveys depicted in Table 1 (from the '80s well into the 2010s). Besides, variations can be limited to only part of the source spectra. Therefore, we decided to estimate the change in the photometric results induced by large variations (worst case scenario) of this source flux, and increase the uncertainties accordingly.

We ran a set of $10^4$ simulations for this source, defining its amplitude at 6.6 GHz from a random Gaussian realization assuming a $\sigma$ equal to 50 per cent of the measured flux in the Fatigoni et al. (2021) catalogue. We assumed this variability level after the analysis ran in Section 4.2 of Harper et al. (2023), where the variability of this source was found to be between 30 and 50 per cent from yearly Planck 28.4 GHz maps. We used the 5C+3.50 broken power spectral model from Harper et al. (2023) to build an estimate for this source intensity at each map frequency. This estimate is then projected onto a map for each frequency assuming the SRT FWHM beam. These maps are smoothed to 1° resolution and transformed to HEALPix maps. These maps, which contain only the signal from 5C+3.50, are then subtracted from the frequency maps. We performed aperture photometry on these final maps, as explained in Section 3.2, and assessed the changes in the extracted flux density for the full M31 aperture. We found that the introduced 50 per cent $\sigma$ variability level translates into a maximum 40 per cent variation in the computed flux density of M31 in the Planck 28.4 GHz band. We increased the flux density error at each frequency by adding in quadrature the uncertainty associated with the variability of this source. Therefore, equation (4) transforms into:

$$ \sigma_{S_v} = \sqrt{\sigma_{\text{AP}}^2 + \text{cal}^2 \cdot S_v^2 + (\sigma_v^{\text{var}})^2 \cdot S_v^2}, \qquad (7) $$

where $\sigma_v^{\text{var}}$ is the variation level found for each band. It is only larger than 10 per cent for WMAP and Planck bands between 20 and 50 GHz, while being negligible for low ($\nu < 10$ GHz) and high ($\nu > 90$ GHz) frequency bands.

## 3.5 Subtraction of CMB anisotropies

The contribution of CMB anisotropies to the global SED of M31 is not negligible, as can be seen in the Planck 143 GHz maps in the bottom row of Fig. 2. In fact, given the relatively large angular size of M31 and its low flux density, CMB anisotropies make an important contribution to the uncertainty in the measured flux densities. As mentioned in Section 2, we used 1° smoothed maps in this study: at that scale, CMB anisotropies show an RMS of $\approx 80$ µK. When translated into flux density, that uncertainty value dominates over the dispersion of the rest of the emission from the background region. It is thus preferable first to subtract the CMB anisotropies and then fit for the other components at the SED level. Several maps tracing CMB anisotropies exist, but we focused on the four official Planck collaboration CMB maps from SMICA, NILC, SEVEM, and COMMANDER methods (Planck Collaboration et al. 2020b).

In order to study to what degree the results presented in this paper could be affected by the selection of the CMB component-separated map, we performed aperture photometry on all CMB maps. First, we used SMICA, following the methodology of B19. We found that the SMICA CMB integrated contribution in the M31 region was the most important of the four maps (+1.0 µK). This estimate (and those for the other maps quoted later on) was computed as the mean temperature within the M31 aperture minus the median of the background region, using the apertures from Section 3.2. COMMANDER is the opposite case: its estimate is the least important of the four (-2.4 µK). NILC and SEVEM yielded similar CMB mean temperatures (-1.3 µK and -0.8 µK), which sit between the SMICA and COMMANDER values, as good trade-off values. We decided to use NILC after performing a visual inspection of the maps (there are visible residuals in SEVEM on the M31 outer ring). Using any of the extreme values (SMICA or COMMANDER) biases the AME estimation to high or low values. This is because the amplitude from the CMB maps within the M31 aperture is comparable to the one from the emission coming from the galaxy itself. This is further discussed in Section 5.1. Finally, in order to
reflect uncertainties associated with the subtraction of the CMB anisotropies, the error bar derived from equation (3) was increased by adding in quadrature the standard deviation of the integrated CMB anisotropy estimates from the previous maps (1.23 µK). In this way, the error budget of the flux density estimates increases mostly in the spectral range where the CMB contribution is important.

## 4 RESULTS

We recover the AME peak flux density parameter, $A_{AME}$, in M31 with a $3.2\sigma$ significance, as can be seen in Table 4 and Fig. 4, $A_{AME} = 1.03 \pm 0.32$ Jy. At 25 GHz, we estimate the AME flux density to be $0.83 \pm 0.32$ Jy, $1.7\sigma$ below the value of $1.45 \pm 0.18$ Jy quoted by B19. The reason for this difference could be that our measured flux densities between 10 and 60 GHz (where AME is expected to be most important) are lower than those from B19, probably due to differences in the filtering of the ancillary data together with the different choices for the CMB map (B19 used SMICA while we are using NILC, as further discussed in Sections 3.5 and 5.1). For example, our photometry estimate for the Planck 28.4 GHz band is 74 per cent than that from B19; correcting for that factor, our $S_{AME}$ estimates at 25 GHz would increase to 1.1 Jy, closer to that from B19. The corner plot showing the parameter posteriors used to build the previous SED is shown in Fig. 5.

Our fitted synchrotron spectral index, $\alpha_{syn} = -0.97^{+0.18}_{-0.23}$ is consistent with those obtained by Berkhuijsen, Beck & Hoernes (2003), $-1.0 \pm 0.2$, and Planck Collaboration et al. (2015b), $-0.92 \pm 0.16$. It is also consistent with the value obtained by B19, $\alpha_{syn} = -1.1 \pm 0.1$, owing to large uncertainty in this parameter due to unresolved parameter degeneracies. This is because of the lack of data with good calibration uncertainties below 10 GHz. It is worth noting that, except for the thermal dust spectral index, the differences between our best-fit parameters and those of B19 are always within the $1\sigma$ level. The study of B19 benefitted from the addition of the SRT data point at 6.6 GHz, which is critical to alleviating the degeneracy between the synchrotron and free-free components. The value of our fit at 6.6 GHz is not too different from the one quoted in the analysis of B19: our flux density estimate value is $1.68^{+0.36}_{-0.34}$ Jy compared to $1.199 \pm 0.087$ Jy from B19. The two are only $1.3\sigma$ away.

Degeneracy with synchrotron emission prevents a precise determination of the level of free-free emission from our component-separation analysis. Actually, Fig. 5 clearly shows that the 1D posterior of EM is consistent with a non-detection of free-free emission, and an upper limit would probably be more appropriate. The 68 per cent upper limit on the free-free flux density at 1 GHz derived from this posterior is 0.52 Jy. This is consistent with the expected value derived from the $H\alpha$ emission. By integration on the map of Finkbeiner (2003) with the aperture and background annulus described in Section 3.2, and using the formalism explained in Dickinson, Davies & Davis (2003), we derive an expected free-free flux density at 1 GHz of $S_{1 GHz}^{ff} = 0.15 \pm 0.02$ Jy. Due to the absorption of $H\alpha$ emission, this value might be regarded as a lower limit. Absorption correction suffers from important uncertainties, related, for instance, with the way absorption measured at other wavelengths is translated into $H\alpha$ absorption. Using the formalism described in Dickinson, Davies & Davis (2003) and the reddening $E(B-V)$ map derived by Planck Collaboration et al. (2014b), we get $S_{1 GHz}^{ff} = 0.19$ Jy. If we instead use the reddening $A_V$ map, also derived by Planck Collaboration et al. (2016a), we get $S_{1 GHz}^{ff} = 0.27$ Jy. These values are both consistent with our derived upper limit.

Another common way of estimating free-free emission relies on measurements of the star formation rate (SFR). Ford et al.

(2013) measured $0.25 \pm 0.05$ M$_\odot$ yr$^{-1}$ from UV and 24-micron data. Equation (6) of Murphy et al. (2012) provides a relation between SFR and spectral luminosity. Using this equation and assuming an average electron temperature $T_e = 8000$ K (typical value from our Galaxy), we get $L_v = 4.92 \times 10^{19}$ W Hz$^{-1}$ at 1 GHz. From this and using a distance to M31 of $d = 785$ pc (McConnachie et al. 2005), we get a flux density at 1 GHz of $S_{1 GHz}^{ff} = 0.67$ Jy, a value that is roughly consistent with our derived upper limit (we have intentionally omitted the statistical error in this estimate as systematic errors associated with the scaling relations used are expected to be larger).

Thus, given that the estimate for the amplitude of free-free is consistent with zero, we repeated our analysis without taking this component into account. The absence of a free-free contribution is compensated by a slightly wider AME component: $W_{AME} = 0.66 \pm 0.18$ now, as compared to $0.58 \pm 0.16$ when free-free was considered.<fnref n="9" /> However, most AME theory models (SPDUST, Ali-Haïmoud, Hirata & Dickinson 2009; Silsbee, Ali-Haïmoud & Hirata 2011) expect values for $W_{AME}$ smaller than 0.6. Previous studies have resolved this issue by adding a second AME component (e.g. Planck Collaboration et al. 2016b). However, in the case of the integrated emission from a galaxy like ours, we expect a large number of AME emitting regions with different spectral parameters within our aperture, leading to a broadening of the integrated AME distribution. We tried to fit a second AME component, but found that the MCMC chains stopped converging and the parameters defining this new component were not properly constrained. In conclusion, we decided to use as the default model in this analysis the one with free-free emission, even though its $\chi_{red}^2$ value is slightly higher than the value when no free-free is considered (0.39 versus 0.61, $\chi_{red}^2$ being systematically smaller than 1 most probably owing to an overestimation of the photometry uncertainties). This decision is based on the fact that we know, from different kinds of observational evidence (e.g. the presence of star formation or $H\alpha$ emission), that free-free emission originating in M31 must exist.

## 4.1 Statistical evaluation of the presence of AME in M31

We repeated the analysis without the AME component to ensure that AME in M31 could not be mistaken for free-free emission. The $\chi_{red}^2$ value obtained in that case is at least a factor 2.5 worse than any of the two other scenarios, where an AME component is introduced (Table 4). $\chi_{red}^2$ values are computed using the median values from the distributions of the parameter posteriors shown in Fig. 5. These median values are those quoted in Table 4. The worst fit is the one for which no AME component is considered.

We used the Akaike information criterion (AIC; Akaike 1974) and the Bayesian information criterion (BIC; Schwarz 1978) to compare the goodness of the fits with and without AME. This was done in order to account for the non-Gaussianity of the likelihood, mostly due to the presence of non-linear parameters in the SED fitting (e.g. the synchrotron spectral index). AIC and BIC are computed as follows:

$$ 
\begin{align} 
\text{AIC} &= 2k - 2 \log \mathcal{L} \\ 
\text{BIC} &= k \log(n) - 2 \log \mathcal{L}, 
\end{align} 
\quad (8) 
$$

where $k$ stands for the number of fitted parameters, while $n$ is the number of data points introduced in the fit. $\log \mathcal{L}$ is the maximum

<footnote n="9">⁹In fact, we can see in the corner plot (Fig. 5) that the free-free amplitude and $W_{AME}$ are slightly anticorrelated. Large values of $W_{AME}$ imply wide AME distributions that could resemble power-law distributions similar to free-free between 10 and 60 GHz, where AME is most important.</footnote>
reflect uncertainties associated with the subtraction of the CMB anisotropies, the error bar derived from equation (3) was increased by adding in quadrature the standard deviation of the integrated CMB anisotropy estimates from the previous maps (1.23 µK). In this way, the error budget of the flux density estimates increases mostly in the spectral range where the CMB contribution is important.

## 4 RESULTS

We recover the AME peak flux density parameter, $A_{AME}$, in M31 with a $3.2\sigma$ significance, as can be seen in Table 4 and Fig. 4, $A_{AME} = 1.03 \pm 0.32$ Jy. At 25 GHz, we estimate the AME flux density to be $0.83 \pm 0.32$ Jy, $1.7\sigma$ below the value of $1.45 \pm 0.18$ Jy quoted by B19. The reason for this difference could be that our measured flux densities between 10 and 60 GHz (where AME is expected to be most important) are lower than those from B19, probably due to differences in the filtering of the ancillary data together with the different choices for the CMB map (B19 used SMICA while we are using NILC, as further discussed in Sections 3.5 and 5.1). For example, our photometry estimate for the Planck 28.4 GHz band is 74 per cent than that from B19; correcting for that factor, our $S_{AME}$ estimates at 25 GHz would increase to 1.1 Jy, closer to that from B19. The corner plot showing the parameter posteriors used to build the previous SED is shown in Fig. 5.

Our fitted synchrotron spectral index, $\alpha_{syn} = -0.97^{+0.18}_{-0.23}$ is consistent with those obtained by Berkhuijsen, Beck & Hoernes (2003), $-1.0 \pm 0.2$, and Planck Collaboration et al. (2015b), $-0.92 \pm 0.16$. It is also consistent with the value obtained by B19, $\alpha_{syn} = -1.1 \pm 0.1$, owing to large uncertainty in this parameter due to unresolved parameter degeneracies. This is because of the lack of data with good calibration uncertainties below 10 GHz. It is worth noting that, except for the thermal dust spectral index, the differences between our best-fit parameters and those of B19 are always within the $1\sigma$ level. The study of B19 benefitted from the addition of the SRT data point at 6.6 GHz, which is critical to alleviating the degeneracy between the synchrotron and free-free components. The value of our fit at 6.6 GHz is not too different from the one quoted in the analysis of B19: our flux density estimate value is $1.68^{+0.36}_{-0.34}$ Jy compared to $1.199 \pm 0.087$ Jy from B19. The two are only $1.3\sigma$ away.

Degeneracy with synchrotron emission prevents a precise determination of the level of free-free emission from our component-separation analysis. Actually, Fig. 5 clearly shows that the 1D posterior of EM is consistent with a non-detection of free-free emission, and an upper limit would probably be more appropriate. The 68 per cent upper limit on the free-free flux density at 1 GHz derived from this posterior is 0.52 Jy. This is consistent with the expected value derived from the $H\alpha$ emission. By integration on the map of Finkbeiner (2003) with the aperture and background annulus described in Section 3.2, and using the formalism explained in Dickinson, Davies & Davis (2003), we derive an expected free-free flux density at 1 GHz of $S_{1 GHz}^{ff} = 0.15 \pm 0.02$ Jy. Due to the absorption of $H\alpha$ emission, this value might be regarded as a lower limit. Absorption correction suffers from important uncertainties, related, for instance, with the way absorption measured at other wavelengths is translated into $H\alpha$ absorption. Using the formalism described in Dickinson, Davies & Davis (2003) and the reddening $E(B-V)$ map derived by Planck Collaboration et al. (2014b), we get $S_{1 GHz}^{ff} = 0.19$ Jy. If we instead use the reddening $A_V$ map, also derived by Planck Collaboration et al. (2016a), we get $S_{1 GHz}^{ff} = 0.27$ Jy. These values are both consistent with our derived upper limit.

Another common way of estimating free-free emission relies on measurements of the star formation rate (SFR). Ford et al.

(2013) measured $0.25 \pm 0.05$ M$_\odot$ yr$^{-1}$ from UV and 24-micron data. Equation (6) of Murphy et al. (2012) provides a relation between SFR and spectral luminosity. Using this equation and assuming an average electron temperature $T_e = 8000$ K (typical value from our Galaxy), we get $L_v = 4.92 \times 10^{19}$ W Hz$^{-1}$ at 1 GHz. From this and using a distance to M31 of $d = 785$ pc (McConnachie et al. 2005), we get a flux density at 1 GHz of $S_{1 GHz}^{ff} = 0.67$ Jy, a value that is roughly consistent with our derived upper limit (we have intentionally omitted the statistical error in this estimate as systematic errors associated with the scaling relations used are expected to be larger).

Thus, given that the estimate for the amplitude of free-free is consistent with zero, we repeated our analysis without taking this component into account. The absence of a free-free contribution is compensated by a slightly wider AME component: $W_{AME} = 0.66 \pm 0.18$ now, as compared to $0.58 \pm 0.16$ when free-free was considered.<fnref n="9" /> However, most AME theory models (SPDUST, Ali-Haïmoud, Hirata & Dickinson 2009; Silsbee, Ali-Haïmoud & Hirata 2011) expect values for $W_{AME}$ smaller than 0.6. Previous studies have resolved this issue by adding a second AME component (e.g. Planck Collaboration et al. 2016b). However, in the case of the integrated emission from a galaxy like ours, we expect a large number of AME emitting regions with different spectral parameters within our aperture, leading to a broadening of the integrated AME distribution. We tried to fit a second AME component, but found that the MCMC chains stopped converging and the parameters defining this new component were not properly constrained. In conclusion, we decided to use as the default model in this analysis the one with free-free emission, even though its $\chi_{red}^2$ value is slightly higher than the value when no free-free is considered (0.39 versus 0.61, $\chi_{red}^2$ being systematically smaller than 1 most probably owing to an overestimation of the photometry uncertainties). This decision is based on the fact that we know, from different kinds of observational evidence (e.g. the presence of star formation or $H\alpha$ emission), that free-free emission originating in M31 must exist.

## 4.1 Statistical evaluation of the presence of AME in M31

We repeated the analysis without the AME component to ensure that AME in M31 could not be mistaken for free-free emission. The $\chi_{red}^2$ value obtained in that case is at least a factor 2.5 worse than any of the two other scenarios, where an AME component is introduced (Table 4). $\chi_{red}^2$ values are computed using the median values from the distributions of the parameter posteriors shown in Fig. 5. These median values are those quoted in Table 4. The worst fit is the one for which no AME component is considered.

We used the Akaike information criterion (AIC; Akaike 1974) and the Bayesian information criterion (BIC; Schwarz 1978) to compare the goodness of the fits with and without AME. This was done in order to account for the non-Gaussianity of the likelihood, mostly due to the presence of non-linear parameters in the SED fitting (e.g. the synchrotron spectral index). AIC and BIC are computed as follows:

$$ 
\begin{align} 
\text{AIC} &= 2k - 2 \log \mathcal{L} \\ 
\text{BIC} &= k \log(n) - 2 \log \mathcal{L}, 
\end{align} 
\quad (8) 
$$

where $k$ stands for the number of fitted parameters, while $n$ is the number of data points introduced in the fit. $\log \mathcal{L}$ is the maximum

<footnote n="9">⁹In fact, we can see in the corner plot (Fig. 5) that the free-free amplitude and $W_{AME}$ are slightly anticorrelated. Large values of $W_{AME}$ imply wide AME distributions that could resemble power-law distributions similar to free-free between 10 and 60 GHz, where AME is most important.</footnote>
<table><caption>**Table 4.** Comparison between the results for the different cases considered, shown in Figs 4 and 6. These values are obtained as the median values from the marginalized 1D distributions, as shown in Fig. 5.</caption><thead><tr><th>Parameter</th><th>Case 1 (all components)</th><th>Case 2 (no AME)</th><th>Case 3 (no free-free)</th><th>From Battistelli et al. (2019)</th></tr></thead><tbody><tr><td>S<sub>1 GHz</sub> (Jy)</td><td>5.91 ± 1.01</td><td>6.90 ± 0.93</td><td>6.50 ± 0.97</td><td>6.97 ± 0.54</td></tr><tr><td>α<sub>syn</sub></td><td>-0.97 ± 0.21</td><td>-0.65 ± 0.08</td><td>-0.89 ± 0.19</td><td>-1.1 ± 0.1</td></tr><tr><td>EM (pc cm<sup>-6</sup>)</td><td>2.72 ± 1.64</td><td>3.53 ± 2.31</td><td>-</td><td>-</td></tr><tr><td>A<sub>AME</sub> (Jy)</td><td>1.03 ± 0.32</td><td>-</td><td>1.20 ± 0.35</td><td>-</td></tr><tr><td>ν<sub>AME</sub> (GHz)</td><td>17.15 ± 3.17</td><td>-</td><td>17.65 ± 3.30</td><td>-</td></tr><tr><td>W<sub>AME</sub></td><td>0.58 ± 0.16</td><td>-</td><td>0.66 ± 0.18</td><td>-</td></tr><tr><td>τ<sub>353</sub></td><td>3.79 ± 0.34</td><td>3.84 ± 0.34</td><td>3.77 ± 0.33</td><td>-</td></tr><tr><td>β<sub>d</sub></td><td>1.71 ± 0.08</td><td>1.79 ± 0.08</td><td>1.68 ± 0.08</td><td>1.49 ± 0.056</td></tr><tr><td>T<sub>d</sub> [K]</td><td>18.49 ± 0.79</td><td>18.02 ± 0.75</td><td>18.69 ± 0.78</td><td>18.8 ± 0.54</td></tr><tr><td>S<sub>1 GHz</sub><sup>ff</sup> (Jy)</td><td>0.39 ± 0.24</td><td>0.51 ± 0.33</td><td>-</td><td>0.33 ± 0.26</td></tr><tr><td>S<sub>AME</sub> @ 25 GHz (Jy)</td><td>0.83 ± 0.32</td><td>-</td><td>1.05 ± 0.35</td><td>1.45 ± 0.18</td></tr><tr><td>S<sub>3000 GHz</sub><sup>dust</sup> (Jy)</td><td>3255 ± 1130</td><td>3161 ± 1098</td><td>3276 ± 1103</td><td>3180 ± 230</td></tr><tr><td>χ<sup>2</sup><sub>red</sub></td><td>0.61</td><td>1.51</td><td>0.39</td><td>0.36</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Top: fitted SED in M31, considering all components (FULL scenario) defined in Section 3.1. Each of the solid lines account for an emission component (blue: synchrotron, green: free-free, red: AME, purple: thermal dust). Grey lines show different realizations of the MCMC. Both fits obtained using the median values from the parameters 1D marginal distributions and the parameter combination returning the maximum likelihood estimate are plotted. Photometry estimates are plotted after applying colour corrections (cc); values for the no-cc case are given in Table 1. The parameters describing the different cases are shown in Table 4. The scaled emission from Berkhuijsen, Beck & Hoernes (2003) shown here is calculated using the spectral index of the full emission, not the one for synchrotron alone. The photometry estimate from B19 at 6.6 GHz is shown for comparison only: it was not used in this fit or in any of the following, although its estimate is consistent with our model estimate at that frequency. Bottom: residuals when comparing the photometry estimates with both fits, units being the errorbars from photometry.</figcaption></figure>
<table><caption>**Table 4.** Comparison between the results for the different cases considered, shown in Figs 4 and 6. These values are obtained as the median values from the marginalized 1D distributions, as shown in Fig. 5.</caption><thead><tr><th>Parameter</th><th>Case 1 (all components)</th><th>Case 2 (no AME)</th><th>Case 3 (no free-free)</th><th>From Battistelli et al. (2019)</th></tr></thead><tbody><tr><td>S<sub>1 GHz</sub> (Jy)</td><td>5.91 ± 1.01</td><td>6.90 ± 0.93</td><td>6.50 ± 0.97</td><td>6.97 ± 0.54</td></tr><tr><td>α<sub>syn</sub></td><td>-0.97 ± 0.21</td><td>-0.65 ± 0.08</td><td>-0.89 ± 0.19</td><td>-1.1 ± 0.1</td></tr><tr><td>EM (pc cm<sup>-6</sup>)</td><td>2.72 ± 1.64</td><td>3.53 ± 2.31</td><td>-</td><td>-</td></tr><tr><td>A<sub>AME</sub> (Jy)</td><td>1.03 ± 0.32</td><td>-</td><td>1.20 ± 0.35</td><td>-</td></tr><tr><td>ν<sub>AME</sub> (GHz)</td><td>17.15 ± 3.17</td><td>-</td><td>17.65 ± 3.30</td><td>-</td></tr><tr><td>W<sub>AME</sub></td><td>0.58 ± 0.16</td><td>-</td><td>0.66 ± 0.18</td><td>-</td></tr><tr><td>τ<sub>353</sub></td><td>3.79 ± 0.34</td><td>3.84 ± 0.34</td><td>3.77 ± 0.33</td><td>-</td></tr><tr><td>β<sub>d</sub></td><td>1.71 ± 0.08</td><td>1.79 ± 0.08</td><td>1.68 ± 0.08</td><td>1.49 ± 0.056</td></tr><tr><td>T<sub>d</sub> [K]</td><td>18.49 ± 0.79</td><td>18.02 ± 0.75</td><td>18.69 ± 0.78</td><td>18.8 ± 0.54</td></tr><tr><td>S<sub>1 GHz</sub><sup>ff</sup> (Jy)</td><td>0.39 ± 0.24</td><td>0.51 ± 0.33</td><td>-</td><td>0.33 ± 0.26</td></tr><tr><td>S<sub>AME</sub> @ 25 GHz (Jy)</td><td>0.83 ± 0.32</td><td>-</td><td>1.05 ± 0.35</td><td>1.45 ± 0.18</td></tr><tr><td>S<sub>3000 GHz</sub><sup>dust</sup> (Jy)</td><td>3255 ± 1130</td><td>3161 ± 1098</td><td>3276 ± 1103</td><td>3180 ± 230</td></tr><tr><td>χ<sup>2</sup><sub>red</sub></td><td>0.61</td><td>1.51</td><td>0.39</td><td>0.36</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top: fitted SED in M31, considering all components (FULL scenario) defined in Section 3.1. Each of the solid lines account for an emission component (blue: synchrotron, green: free-free, red: AME, purple: thermal dust). Grey lines show different realizations of the MCMC. Both fits obtained using the median values from the parameters 1D marginal distributions and the parameter combination returning the maximum likelihood estimate are plotted. Photometry estimates are plotted after applying colour corrections (cc); values for the no-cc case are given in Table 1. The parameters describing the different cases are shown in Table 4. The scaled emission from Berkhuijsen, Beck & Hoernes (2003) shown here is calculated using the spectral index of the full emission, not the one for synchrotron alone. The photometry estimate from B19 at 6.6 GHz is shown for comparison only: it was not used in this fit or in any of the following, although its estimate is consistent with our model estimate at that frequency. Bottom: residuals when comparing the photometry estimates with both fits, units being the errorbars from photometry.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Corner plot example, from the FULL case (i.e. the one also shown in Fig. 4). We can see how the most important degeneracies (aside from the ones between the thermal dust components, which are related to the modelling) are those between the AME amplitude ($A_{AME}$) and the free-free one (EM), but EM is also degenerate with $W_{AME}$. This degeneracy arises because a wide enough AME component would resemble the behaviour of a flat power law (similar to a free-free component) for those frequencies where free-free is most important (between 5 and 100 GHz). Solid line contours encompass 1, 2, and $3\sigma$ levels.</figcaption></figure>

log-likelihood value computed within the MCMC, as explained in Section 3.3 and equation (5). Larger differences on the AIC or BIC estimates from two models will imply a heavier preference of the data for the model with the lowest value. The main difference between the two is that the BIC penalizes the addition of parameters more heavily than the AIC. We compare the case with all the components and the case without AME, as we previously explained that the case with no free-free component is not physical:

$$AIC_{\text{all}} = 23.6; \quad AIC_{\text{no AME}} = 31.9;$$

$$BIC_{\text{all}} = 32.5; \quad BIC_{\text{no AME}} = 37.9.$$

According to the AIC, the relative likelihood of the model without an AME component compared to the model with all components considered is:

$$\exp((AIC_{\text{all}} - AIC_{\text{no AME}})/2) = 0.0016,$$

and the model with AME is strongly preferred. When using the BIC, we found the relative likelihood of the model without AME to be:

$$\exp((BIC_{\text{all}} - BIC_{\text{no AME}})/2) = 0.0067,$$

over a factor 4 larger. This is consistent with BIC penalizing the larger number of parameters considered in the case with all components, compared to that without AME. However, the model without AME
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Corner plot example, from the FULL case (i.e. the one also shown in Fig. 4). We can see how the most important degeneracies (aside from the ones between the thermal dust components, which are related to the modelling) are those between the AME amplitude ($A_{AME}$) and the free-free one (EM), but EM is also degenerate with $W_{AME}$. This degeneracy arises because a wide enough AME component would resemble the behaviour of a flat power law (similar to a free-free component) for those frequencies where free-free is most important (between 5 and 100 GHz). Solid line contours encompass 1, 2, and $3\sigma$ levels.</figcaption></figure>

log-likelihood value computed within the MCMC, as explained in Section 3.3 and equation (5). Larger differences on the AIC or BIC estimates from two models will imply a heavier preference of the data for the model with the lowest value. The main difference between the two is that the BIC penalizes the addition of parameters more heavily than the AIC. We compare the case with all the components and the case without AME, as we previously explained that the case with no free-free component is not physical:

$$AIC_{\text{all}} = 23.6; \quad AIC_{\text{no AME}} = 31.9;$$

$$BIC_{\text{all}} = 32.5; \quad BIC_{\text{no AME}} = 37.9.$$

According to the AIC, the relative likelihood of the model without an AME component compared to the model with all components considered is:

$$\exp((AIC_{\text{all}} - AIC_{\text{no AME}})/2) = 0.0016,$$

and the model with AME is strongly preferred. When using the BIC, we found the relative likelihood of the model without AME to be:

$$\exp((BIC_{\text{all}} - BIC_{\text{no AME}})/2) = 0.0067,$$

over a factor 4 larger. This is consistent with BIC penalizing the larger number of parameters considered in the case with all components, compared to that without AME. However, the model without AME
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Same M31 SED as in Fig. 4, but now without a free-free component (top panel, noFF), and finally without an AME component (bottom panel, noAME). The parameters describing the different cases are shown in Table 4.</figcaption></figure>

is still below 0.01 times as probable as that with AME, pointing to a strong preference for the latter.

## 4.2 The global properties of M31 compared to the Milky Way

As pointed out in Section 1, M31 is the galaxy most similar to the MW in the Local Group. Therefore, it is straightforward to compare the properties between the two, especially for AME, which is the main focus of our study. First, we computed the AME emissivity for M31 as done in Fernández-Torreiro et al. (2023), i.e. as the ratio between the AME intensity at 28.4 GHz (in temperature units) and the dust intensity at 100 µm from the fit:

$$ \epsilon_{\text{AME}}^{28.4 \text{ GHz}} = \frac{T_{\text{AME}}^{28.4 \text{ GHz}}}{I_{100 \text{ µm}}} = \frac{\frac{c^2}{2k_B V^2} S_{\text{AME}}^{28.4 \text{ GHz}}}{S_{100 \text{ µm}}}. \quad (9) $$

This ratio is intended to cancel the dependence on the column density present in both AME and thermal dust emission, thus allowing us to compare regions with extremely different morphologies. As pointed out by Tibbs, Paladini & Dickinson (2012), this cancellation may not be perfect as this calculation of the AME emissivity is sensitive to the dust temperature. However, thanks to the similarity between dust temperatures in M31 and in the MW, this should not be a problem in our comparison.

Differences between B19 and this work on dust parameters arise because of the different surveys used (Herschel and COBE-DIRBE, respectively). We find $\epsilon_{\text{AME}}^{28.4 \text{ GHz}} = 9.6 \pm 3.1 \text{ µK MJy}^{-1}$ sr and compare it to the results for the MW both from the Planck

collaboration (Planck Collaboration et al. 2015a) and also Section 4.2 of Fernández-Torreiro et al. (2023). Planck Collaboration et al. (2015a) quoted a $9.8 \pm 0.5 \text{ µK MJy}^{-1}$ sr value<fnref n="10" /> for the MW, fully consistent with our M31 result, while Fernández-Torreiro et al. (2023) returns $\epsilon_{\text{AME}}^{28.4 \text{ GHz}} = 8.8 \pm 3.8 \text{ µK MJy}^{-1}$ sr, which is a lower value, but still consistent because of the large uncertainties.<fnref n="11" /> Both values are lower than those previously found in several other works (see e.g. Table 3 of Dickinson et al. 2018), but the AME emissivity can show large variations (Davies et al. 2006 found variations as large as a factor 2, for example). Those variations would probably be embedded within the integrated spectrum of the galaxy, effectively decreasing the AME emissivity value with respect to that found in compact sources.

When computing the AME fraction at 28.4 GHz, we find $S_{\text{AME}}^{28.4 \text{ GHz}}/S_{\text{total}}^{28.4 \text{ GHz}} = 0.54 \pm 0.17$. This is slightly greater than the values obtained by Fernández-Torreiro et al. (2023) and Planck Collaboration et al. (2015a), $0.46 \pm 0.08$ and $0.45 \pm 0.01$, respectively, although it is consistent to within $1\sigma$, owing to the greater uncertainty in the M31 estimate. Finally, when focusing on the synchrotron spectral index, we find consistent estimates between M31 and the MW of $-0.97^{+0.18}_{-0.23}$ and $-0.94 \pm 0.10$ (Fernández-Torreiro et al. 2023), respectively.

## 4.3 QUIJOTE-MFI M31 upper limits on polarization

As explained in Section 2.1, the M31 field has the best sensitivity values (shown in Table 2) of all the fields studied by QUIJOTE-MFI so far, being slightly better than those from the Haze (Guidi et al. 2023) and W49, W51, and IC443 (Tramonte et al. 2023) fields. This is because of the long integration times invested in those regions. Stokes Q and U maps of the M31 field are displayed in Fig. 7.

Given the lack of any detectable polarized signal towards M31 on the QUIJOTE-MFI maps, we derived upper limits on its global polarized flux density. Using the same photometric approach and apertures as for intensity (Section 3.2), the measured Stokes Q and U parameters are $-0.046 \pm 0.061$ and $0.115 \pm 0.067$ Jy at 11 GHz, and $-0.122 \pm 0.050$ and $0.088 \pm 0.061$ Jy at 13 GHz. Except for the Q value at 13 GHz showing a marginal detection, they are statistically consistent with zero. Taking the measured values as upper limits, the noise unbiased<fnref n="12" /> polarized intensity estimates, P, are $<0.22$ and $<0.18$ Jy, respectively at 11 and 13 GHz for a 95 per cent confidence limit (C.L.). The flux density estimates at 11 and 13 GHz are $1.79 \pm 0.30$ and $1.94 \pm 0.36$ Jy (after subtracting both the point sources and the CMB, as shown in Table 1), so the polarization fractions for the integrated M31 flux density are below 12.5 and 9.6 per cent, respectively, for a 95 per cent C.L. Considering the AME contribution to the total flux in intensity, the polarization fraction upper limits derived from this analysis are $\le 30$ per cent at 11 GHz and $\le 20$ per cent at 13 GHz. Although these upper limits are an order of magnitude greater than the most stringent ones obtained from the analysis of AME sources in the MW (e.g. Génova-Santos et al. 2017; Tramonte et al. 2023; González-González et al., in preparation), this is the first time that a limit on the AME polarization has been obtained for an extragalactic object. It must be noted, though, that a potential underlying AME polarization in M31 in

<footnote n="10">¹⁰There seems to be a typo on this Planck paper, as in Section 7.3 the units for the emissivity are quoted as $\text{mK MJy}^{-1}$ instead of $\text{µK MJy}^{-1}$ sr.</footnote>

<footnote n="11">¹¹This estimate increases when only those pixels with strong AME detections are used, up to $\epsilon_{\text{AME}}^{28.4 \text{ GHz}} = 11.6 \pm 3.5 \text{ µK MJy}^{-1}$ sr.</footnote>

<footnote n="12">¹²Computed following Vaillancourt (2006) and Rubiño-Martín et al. (2012).</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same M31 SED as in Fig. 4, but now without a free-free component (top panel, noFF), and finally without an AME component (bottom panel, noAME). The parameters describing the different cases are shown in Table 4.</figcaption></figure>

is still below 0.01 times as probable as that with AME, pointing to a strong preference for the latter.

## 4.2 The global properties of M31 compared to the Milky Way

As pointed out in Section 1, M31 is the galaxy most similar to the MW in the Local Group. Therefore, it is straightforward to compare the properties between the two, especially for AME, which is the main focus of our study. First, we computed the AME emissivity for M31 as done in Fernández-Torreiro et al. (2023), i.e. as the ratio between the AME intensity at 28.4 GHz (in temperature units) and the dust intensity at 100 µm from the fit:

$$ \epsilon_{\text{AME}}^{28.4 \text{ GHz}} = \frac{T_{\text{AME}}^{28.4 \text{ GHz}}}{I_{100 \text{ µm}}} = \frac{\frac{c^2}{2k_B V^2} S_{\text{AME}}^{28.4 \text{ GHz}}}{S_{100 \text{ µm}}}. \quad (9) $$

This ratio is intended to cancel the dependence on the column density present in both AME and thermal dust emission, thus allowing us to compare regions with extremely different morphologies. As pointed out by Tibbs, Paladini & Dickinson (2012), this cancellation may not be perfect as this calculation of the AME emissivity is sensitive to the dust temperature. However, thanks to the similarity between dust temperatures in M31 and in the MW, this should not be a problem in our comparison.

Differences between B19 and this work on dust parameters arise because of the different surveys used (Herschel and COBE-DIRBE, respectively). We find $\epsilon_{\text{AME}}^{28.4 \text{ GHz}} = 9.6 \pm 3.1 \text{ µK MJy}^{-1}$ sr and compare it to the results for the MW both from the Planck

collaboration (Planck Collaboration et al. 2015a) and also Section 4.2 of Fernández-Torreiro et al. (2023). Planck Collaboration et al. (2015a) quoted a $9.8 \pm 0.5 \text{ µK MJy}^{-1}$ sr value<fnref n="10" /> for the MW, fully consistent with our M31 result, while Fernández-Torreiro et al. (2023) returns $\epsilon_{\text{AME}}^{28.4 \text{ GHz}} = 8.8 \pm 3.8 \text{ µK MJy}^{-1}$ sr, which is a lower value, but still consistent because of the large uncertainties.<fnref n="11" /> Both values are lower than those previously found in several other works (see e.g. Table 3 of Dickinson et al. 2018), but the AME emissivity can show large variations (Davies et al. 2006 found variations as large as a factor 2, for example). Those variations would probably be embedded within the integrated spectrum of the galaxy, effectively decreasing the AME emissivity value with respect to that found in compact sources.

When computing the AME fraction at 28.4 GHz, we find $S_{\text{AME}}^{28.4 \text{ GHz}}/S_{\text{total}}^{28.4 \text{ GHz}} = 0.54 \pm 0.17$. This is slightly greater than the values obtained by Fernández-Torreiro et al. (2023) and Planck Collaboration et al. (2015a), $0.46 \pm 0.08$ and $0.45 \pm 0.01$, respectively, although it is consistent to within $1\sigma$, owing to the greater uncertainty in the M31 estimate. Finally, when focusing on the synchrotron spectral index, we find consistent estimates between M31 and the MW of $-0.97^{+0.18}_{-0.23}$ and $-0.94 \pm 0.10$ (Fernández-Torreiro et al. 2023), respectively.

## 4.3 QUIJOTE-MFI M31 upper limits on polarization

As explained in Section 2.1, the M31 field has the best sensitivity values (shown in Table 2) of all the fields studied by QUIJOTE-MFI so far, being slightly better than those from the Haze (Guidi et al. 2023) and W49, W51, and IC443 (Tramonte et al. 2023) fields. This is because of the long integration times invested in those regions. Stokes Q and U maps of the M31 field are displayed in Fig. 7.

Given the lack of any detectable polarized signal towards M31 on the QUIJOTE-MFI maps, we derived upper limits on its global polarized flux density. Using the same photometric approach and apertures as for intensity (Section 3.2), the measured Stokes Q and U parameters are $-0.046 \pm 0.061$ and $0.115 \pm 0.067$ Jy at 11 GHz, and $-0.122 \pm 0.050$ and $0.088 \pm 0.061$ Jy at 13 GHz. Except for the Q value at 13 GHz showing a marginal detection, they are statistically consistent with zero. Taking the measured values as upper limits, the noise unbiased<fnref n="12" /> polarized intensity estimates, P, are $<0.22$ and $<0.18$ Jy, respectively at 11 and 13 GHz for a 95 per cent confidence limit (C.L.). The flux density estimates at 11 and 13 GHz are $1.79 \pm 0.30$ and $1.94 \pm 0.36$ Jy (after subtracting both the point sources and the CMB, as shown in Table 1), so the polarization fractions for the integrated M31 flux density are below 12.5 and 9.6 per cent, respectively, for a 95 per cent C.L. Considering the AME contribution to the total flux in intensity, the polarization fraction upper limits derived from this analysis are $\le 30$ per cent at 11 GHz and $\le 20$ per cent at 13 GHz. Although these upper limits are an order of magnitude greater than the most stringent ones obtained from the analysis of AME sources in the MW (e.g. Génova-Santos et al. 2017; Tramonte et al. 2023; González-González et al., in preparation), this is the first time that a limit on the AME polarization has been obtained for an extragalactic object. It must be noted, though, that a potential underlying AME polarization in M31 in

<footnote n="10">¹⁰There seems to be a typo on this Planck paper, as in Section 7.3 the units for the emissivity are quoted as $\text{mK MJy}^{-1}$ instead of $\text{µK MJy}^{-1}$ sr.</footnote>

<footnote n="11">¹¹This estimate increases when only those pixels with strong AME detections are used, up to $\epsilon_{\text{AME}}^{28.4 \text{ GHz}} = 11.6 \pm 3.5 \text{ µK MJy}^{-1}$ sr.</footnote>

<footnote n="12">¹²Computed following Vaillancourt (2006) and Rubiño-Martín et al. (2012).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 7.** QUIJOTE-MFI maps at 11.2 and 12.9 GHz (top and bottom rows, respectively) for the M31 field in intensity (left column) and Stokes parameters *Q* (centre) and *U* (right). The solid lines show the M31 (ellipse) and background (rectangular) apertures defined in Section 3.2 and used in Section 4. These maps are the raw ones after combining the data from the WS and raster observations: we have not subtracted neither the CMB nor point sources.</figcaption></figure>

small angular scales may be largely smeared out in our aperture due to the mixing of different polarization orientations. This would also be the case for a possible polarization of the synchrotron emission, that in specific regions of M31 is known to reach values of ~60 per cent (Berkhuijsen, Beck & Hoernes 2003). In any case, owing to the sensitivity of our polarization data, and to the uncertainty in the determination of the integrated total-intensity flux density of the synchrotron emission, our derived constraints on the global synchrotron polarization fraction are above this level.

# 5 ROBUSTNESS TESTS

## 5.1 Impact of the selected CMB map on the AME amplitude

As first mentioned in Section 3.5, different CMB maps yield different integrated measurements on the M31 aperture. These affect the photometry flux densities computed to build the SED as in equation (2). We mentioned that the measurement from SMICA is the largest of the four maps, which implies a decrease in M31 flux densities at microwave frequencies ($v \in (40, 90)$ GHz). This produces a steeper spectrum in this frequency range, and an AME component is preferred over free-free or synchrotron. In the opposite scenario, the measurement from COMMANDER is lower, thus flattening the SED at those frequencies and decreasing AME amplitude.

We ran the same analysis as in Section 4 after using SMICA and COMMANDER instead of NILC to subtract the CMB anisotropies. AME amplitude increases to $A_{\text{AME}} = 1.17 \pm 0.29$ Jy when the SMICA map is used, and decreases to $A_{\text{AME}} = 0.92 \pm 0.34$ Jy when using COMMANDER. Thus, the AME significance changes from $4.0\sigma$ to $2.7\sigma$ depending on the choice of the CMB map. Besides, the higher $A_{\text{AME}}$ value when subtracting the signal from SMICA implies a lower free-free amplitude and steeper synchrotron index. In the COMMANDER case, the situation is the opposite: the $A_{\text{AME}}$ decrease implies a larger free-free component and a flatter synchrotron

spectrum. Because of the differences between the two cases, we decided to use the NILC CMB map, which returns a $3.2\sigma$ $A_{\text{AME}}$ significance instead, as a trade-off scenario.

## 5.2 QUIJOTE-MFI role in the fitting

We repeated the former analyses without taking into account the data from QUIJOTE-MFI. In this way, we can assess its constraining power. We show the results in Table 5 and the corresponding SEDs for the FULL (all emission components considered) case in Fig. 8. The significance of the measurement of $A_{\text{AME}}$ is slightly reduced to $3.1\sigma$ when the QUIJOTE-MFI data are not used: it is now $1.03 \pm 0.33$ Jy compared to the previous value of $1.03 \pm 0.32$ Jy. However, we find that the $v_{\text{AME}}$ posterior is now not constrained at lower values. This implies an increase on its uncertainty, so now $v_{\text{AME}} = 18.3 \pm 5.8$ GHz, compared to $17.2 \pm 3.2$ GHz when using QUIJOTE-MFI data. This is directly related to the absence of data in that frequency range (5–20 GHz). The larger uncertainty also implies larger differences between the preferred models when using the max-likelihood solution or that from the 1D distributions median values. The change is not as important for the AME width, $W_{\text{AME}}$, as it is $0.54 \pm 0.20$ now compared to the previous value of $0.58 \pm 0.16$. Therefore, QUIJOTE-MFI data addition significantly reduces (almost by a factor 2) the uncertainty of the $v_{\text{AME}}$ parameter. Finally, we can see that the model providing the worst fit to the data is again the one trying to account for all emissions with just synchrotron and free-free components and no AME component.

## 5.3 Source subtraction effect

We repeated our analysis of the maps with no point source subtraction applied, but keeping QUIJOTE-MFI data. This was done in order to assess the impact of this subtraction on the level of AME. Results for this case are shown in Table 6. The FULL
<figure><img src="image_9.png" /><figcaption>**Figure 7.** QUIJOTE-MFI maps at 11.2 and 12.9 GHz (top and bottom rows, respectively) for the M31 field in intensity (left column) and Stokes parameters *Q* (centre) and *U* (right). The solid lines show the M31 (ellipse) and background (rectangular) apertures defined in Section 3.2 and used in Section 4. These maps are the raw ones after combining the data from the WS and raster observations: we have not subtracted neither the CMB nor point sources.</figcaption></figure>

small angular scales may be largely smeared out in our aperture due to the mixing of different polarization orientations. This would also be the case for a possible polarization of the synchrotron emission, that in specific regions of M31 is known to reach values of ~60 per cent (Berkhuijsen, Beck & Hoernes 2003). In any case, owing to the sensitivity of our polarization data, and to the uncertainty in the determination of the integrated total-intensity flux density of the synchrotron emission, our derived constraints on the global synchrotron polarization fraction are above this level.

# 5 ROBUSTNESS TESTS

## 5.1 Impact of the selected CMB map on the AME amplitude

As first mentioned in Section 3.5, different CMB maps yield different integrated measurements on the M31 aperture. These affect the photometry flux densities computed to build the SED as in equation (2). We mentioned that the measurement from SMICA is the largest of the four maps, which implies a decrease in M31 flux densities at microwave frequencies ($v \in (40, 90)$ GHz). This produces a steeper spectrum in this frequency range, and an AME component is preferred over free-free or synchrotron. In the opposite scenario, the measurement from COMMANDER is lower, thus flattening the SED at those frequencies and decreasing AME amplitude.

We ran the same analysis as in Section 4 after using SMICA and COMMANDER instead of NILC to subtract the CMB anisotropies. AME amplitude increases to $A_{\text{AME}} = 1.17 \pm 0.29$ Jy when the SMICA map is used, and decreases to $A_{\text{AME}} = 0.92 \pm 0.34$ Jy when using COMMANDER. Thus, the AME significance changes from $4.0\sigma$ to $2.7\sigma$ depending on the choice of the CMB map. Besides, the higher $A_{\text{AME}}$ value when subtracting the signal from SMICA implies a lower free-free amplitude and steeper synchrotron index. In the COMMANDER case, the situation is the opposite: the $A_{\text{AME}}$ decrease implies a larger free-free component and a flatter synchrotron

spectrum. Because of the differences between the two cases, we decided to use the NILC CMB map, which returns a $3.2\sigma$ $A_{\text{AME}}$ significance instead, as a trade-off scenario.

## 5.2 QUIJOTE-MFI role in the fitting

We repeated the former analyses without taking into account the data from QUIJOTE-MFI. In this way, we can assess its constraining power. We show the results in Table 5 and the corresponding SEDs for the FULL (all emission components considered) case in Fig. 8. The significance of the measurement of $A_{\text{AME}}$ is slightly reduced to $3.1\sigma$ when the QUIJOTE-MFI data are not used: it is now $1.03 \pm 0.33$ Jy compared to the previous value of $1.03 \pm 0.32$ Jy. However, we find that the $v_{\text{AME}}$ posterior is now not constrained at lower values. This implies an increase on its uncertainty, so now $v_{\text{AME}} = 18.3 \pm 5.8$ GHz, compared to $17.2 \pm 3.2$ GHz when using QUIJOTE-MFI data. This is directly related to the absence of data in that frequency range (5–20 GHz). The larger uncertainty also implies larger differences between the preferred models when using the max-likelihood solution or that from the 1D distributions median values. The change is not as important for the AME width, $W_{\text{AME}}$, as it is $0.54 \pm 0.20$ now compared to the previous value of $0.58 \pm 0.16$. Therefore, QUIJOTE-MFI data addition significantly reduces (almost by a factor 2) the uncertainty of the $v_{\text{AME}}$ parameter. Finally, we can see that the model providing the worst fit to the data is again the one trying to account for all emissions with just synchrotron and free-free components and no AME component.

## 5.3 Source subtraction effect

We repeated our analysis of the maps with no point source subtraction applied, but keeping QUIJOTE-MFI data. This was done in order to assess the impact of this subtraction on the level of AME. Results for this case are shown in Table 6. The FULL
<table><caption>**Table 5.** Same as Table 4, but without taking into account QUIJOTE-MFI data.</caption><thead><tr><th>Parameter</th><th>Case 1 (all components)</th><th>Case 2 (no AME)</th><th>Case 3 (no free-free)</th><th>From Battistelli et al. (2019)</th></tr></thead><tbody><tr><td>S<sub>1 GHz</sub> (Jy)</td><td>5.92 ± 1.01</td><td>7.10 ± 0.95</td><td>6.51 ± 0.99</td><td>6.97 ± 0.54</td></tr><tr><td>α<sub>syn</sub></td><td>-0.99 ± 0.21</td><td>-0.65 ± 0.08</td><td>-0.89 ± 0.19</td><td>-1.1 ± 0.1</td></tr><tr><td>EM (pc cm<sup>-6</sup>)</td><td>2.77 ± 1.65</td><td>3.19 ± 2.15</td><td>-</td><td>-</td></tr><tr><td>A<sub>AME</sub> (Jy)</td><td>1.03 ± 0.33</td><td>-</td><td>1.20 ± 0.35</td><td>-</td></tr><tr><td>ν<sub>AME</sub> (GHz)</td><td>18.29 ± 5.84</td><td>-</td><td>18.39 ± 5.79</td><td>-</td></tr><tr><td>W<sub>AME</sub></td><td>0.54 ± 0.20</td><td>-</td><td>0.63 ± 0.22</td><td>-</td></tr><tr><td>τ<sub>353</sub></td><td>3.79 ± 0.34</td><td>3.85 ± 0.35</td><td>3.76 ± 0.33</td><td>-</td></tr><tr><td>β<sub>d</sub></td><td>1.71 ± 0.08</td><td>1.79 ± 0.08</td><td>1.68 ± 0.08</td><td>1.49 ± 0.056</td></tr><tr><td>T<sub>d</sub> [K]</td><td>18.49 ± 0.78</td><td>18.01 ± 0.75</td><td>18.72 ± 0.79</td><td>18.8 ± 0.54</td></tr><tr><td>S<sup>ff</sup><sub>1 GHz</sub> (Jy)</td><td>0.40 ± 0.24</td><td>0.46 ± 0.31</td><td>-</td><td>0.33 ± 0.26</td></tr><tr><td>S<sup>AME</sup> @ 25 GHz (Jy)</td><td>0.87 ± 0.42</td><td>-</td><td>1.06 ± 0.41</td><td>1.45 ± 0.18</td></tr><tr><td>S<sup>dust</sup><sub>3000 GHz</sub> (Jy)</td><td>3219 ± 1104</td><td>3155 ± 1102</td><td>3284 ± 1113</td><td>3180 ± 230</td></tr><tr><td>χ<sup>red</sup></td><td>0.62</td><td>1.53</td><td>0.39</td><td>0.36</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Top panel: same as Fig. 4, but now without taking into account QUIJOTE-MFI data (noMFI). We can see how the dispersion of the model is now much greater between 1 and 20 GHz, as expected. This causes a larger disagreement between the median and maximum likelihood fits than in the rest of the scenarios. Bottom panel: same, but now including QUIJOTE-MFI data and without subtracting point sources emission (noCSS).</figcaption></figure>

case is also shown in the bottom panel of Fig. 8. As it was the case when discarding QUIJOTE-MFI data, the uncertainties on the AME parameters increased, thus reducing the significance of those parameters. There are slight changes in some other parameters. In particular, the synchrotron becomes flatter, and as a result there is more synchrotron at microwave frequencies. The level of free-free is higher: as a consequence of both changes, the AME amplitude is

lower. Therefore, the significance of $A_{AME}$ is now lower too, down to $2.3\sigma$. The $\nu_{AME}$ and $W_{AME}$ values do not change significantly with respect to the main case presented in Section 4: $19.5 \pm 4.3$ GHz and $0.58 \pm 0.26$, respectively, compared to $17.2 \pm 3.2$ GHz and $0.58 \pm 0.16$. The changes of all parameters are always within $1\sigma$.

## 5.4 Impact from changes in the background aperture

Finally, we repeated the analysis using the same elliptical apertures used in Harper et al. (2023), which are also similar to the ones introduced in Planck Collaboration et al. (2015b). These regions have major semiaxes equal to 100, 110, and 154 arcmin, the last two being the inner and outer radii of the ring enclosing the background region. The minor-to-major ratio is 0.7 for all three, with an east-to-north position angle equal to -45°. When comparing with Planck Collaboration et al. (2015b), we are introducing the MFI data and subtracting the CMB using NILC as a template (plus adding an uncertainty term to account for the dispersion in the values of the CMB templates, as explained in Section 3.5). Our main difference with Harper et al. (2023) is the addition of the QUIJOTE-MFI raster data presented in Section 2.1, while we lack the point from C-BASS at 4.76 GHz. That paper also used SMICA instead of NILC to account for CMB anisotropies; we found that SMICA systematically returned higher estimates for the AME than the rest of maps (following the discussion in Section 3.5). We are also not keeping the point source subtraction because of the paucity of objects in Fatigoni et al. (2021) outside the SRT M31 field. Results are shown in Table 7 (where we show the case from Harper et al. (2023) with no subtraction of sources, as it is the most similar one to our analysis) and the corresponding SEDs in Fig. 9. We recover an AME estimate with reduced significance ($2.2\sigma$) only: $0.93 \pm 0.43$ Jy. Note, though, that the differences between our parameters and those of Harper et al. (2023) are always below the $1\sigma$ level.

We find that the result in our study for the synchrotron index, $\alpha_{syn}$, is consistent with those of Berkhuijsen, Beck & Hoernes (2003) and Planck Collaboration et al. (2015b) ($-0.90 \pm 0.20$ versus $-1.0 \pm 0.2$ and $-0.92 \pm 0.16$, respectively), but steeper than that of Harper et al. (2023) ($-0.63 \pm 0.05$). The synchrotron amplitude, $A_{1 GHz}$ is consistent between this study and Planck Collaboration et al. (2015b) and Harper et al. (2023): $8.6 \pm 1.7$ Jy versus $9.5 \pm 1.1$ and $9.7 \pm 0.6$, respectively. Our free-free estimate is larger, although consistent with a non-detection to within $2\sigma$ too: $5.7 \pm 3.4$ pc cm$^{-6}$ versus $1.8 \pm 1.3$ and $1.1 \pm 0.8$ pc cm$^{-6}$. We get a similar consistency in the AME amplitude with Planck Collaboration et al. (2015b)
<table><caption>**Table 5.** Same as Table 4, but without taking into account QUIJOTE-MFI data.</caption><thead><tr><th>Parameter</th><th>Case 1 (all components)</th><th>Case 2 (no AME)</th><th>Case 3 (no free-free)</th><th>From Battistelli et al. (2019)</th></tr></thead><tbody><tr><td>S<sub>1 GHz</sub> (Jy)</td><td>5.92 ± 1.01</td><td>7.10 ± 0.95</td><td>6.51 ± 0.99</td><td>6.97 ± 0.54</td></tr><tr><td>α<sub>syn</sub></td><td>-0.99 ± 0.21</td><td>-0.65 ± 0.08</td><td>-0.89 ± 0.19</td><td>-1.1 ± 0.1</td></tr><tr><td>EM (pc cm<sup>-6</sup>)</td><td>2.77 ± 1.65</td><td>3.19 ± 2.15</td><td>—</td><td>—</td></tr><tr><td>A<sub>AME</sub> (Jy)</td><td>1.03 ± 0.33</td><td>—</td><td>1.20 ± 0.35</td><td>—</td></tr><tr><td>ν<sub>AME</sub> (GHz)</td><td>18.29 ± 5.84</td><td>—</td><td>18.39 ± 5.79</td><td>—</td></tr><tr><td>W<sub>AME</sub></td><td>0.54 ± 0.20</td><td>—</td><td>0.63 ± 0.22</td><td>—</td></tr><tr><td>τ<sub>353</sub></td><td>3.79 ± 0.34</td><td>3.85 ± 0.35</td><td>3.76 ± 0.33</td><td>—</td></tr><tr><td>β<sub>d</sub></td><td>1.71 ± 0.08</td><td>1.79 ± 0.08</td><td>1.68 ± 0.08</td><td>1.49 ± 0.056</td></tr><tr><td>T<sub>d</sub> [K]</td><td>18.49 ± 0.78</td><td>18.01 ± 0.75</td><td>18.72 ± 0.79</td><td>18.8 ± 0.54</td></tr><tr><td>S<sup>ff</sup><sub>1 GHz</sub> (Jy)</td><td>0.40 ± 0.24</td><td>0.46 ± 0.31</td><td>—</td><td>0.33 ± 0.26</td></tr><tr><td>S<sup>AME</sup> @ 25 GHz (Jy)</td><td>0.87 ± 0.42</td><td>—</td><td>1.06 ± 0.41</td><td>1.45 ± 0.18</td></tr><tr><td>S<sup>dust</sup><sub>3000 GHz</sub> (Jy)</td><td>3219 ± 1104</td><td>3155 ± 1102</td><td>3284 ± 1113</td><td>3180 ± 230</td></tr><tr><td>χ<sup>red</sup></td><td>0.62</td><td>1.53</td><td>0.39</td><td>0.36</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Top panel: same as Fig. 4, but now without taking into account QUIJOTE-MFI data (noMFI). We can see how the dispersion of the model is now much greater between 1 and 20 GHz, as expected. This causes a larger disagreement between the median and maximum likelihood fits than in the rest of the scenarios. Bottom panel: same, but now including QUIJOTE-MFI data and without subtracting point sources emission (noCSS).</figcaption></figure>

case is also shown in the bottom panel of Fig. 8. As it was the case when discarding QUIJOTE-MFI data, the uncertainties on the AME parameters increased, thus reducing the significance of those parameters. There are slight changes in some other parameters. In particular, the synchrotron becomes flatter, and as a result there is more synchrotron at microwave frequencies. The level of free-free is higher: as a consequence of both changes, the AME amplitude is

lower. Therefore, the significance of $A_{AME}$ is now lower too, down to $2.3\sigma$. The $\nu_{AME}$ and $W_{AME}$ values do not change significantly with respect to the main case presented in Section 4: $19.5 \pm 4.3$ GHz and $0.58 \pm 0.26$, respectively, compared to $17.2 \pm 3.2$ GHz and $0.58 \pm 0.16$. The changes of all parameters are always within $1\sigma$.

## 5.4 Impact from changes in the background aperture

Finally, we repeated the analysis using the same elliptical apertures used in Harper et al. (2023), which are also similar to the ones introduced in Planck Collaboration et al. (2015b). These regions have major semiaxes equal to 100, 110, and 154 arcmin, the last two being the inner and outer radii of the ring enclosing the background region. The minor-to-major ratio is 0.7 for all three, with an east-to-north position angle equal to $-45^\circ$. When comparing with Planck Collaboration et al. (2015b), we are introducing the MFI data and subtracting the CMB using NILC as a template (plus adding an uncertainty term to account for the dispersion in the values of the CMB templates, as explained in Section 3.5). Our main difference with Harper et al. (2023) is the addition of the QUIJOTE-MFI raster data presented in Section 2.1, while we lack the point from C-BASS at 4.76 GHz. That paper also used SMICA instead of NILC to account for CMB anisotropies; we found that SMICA systematically returned higher estimates for the AME than the rest of maps (following the discussion in Section 3.5). We are also not keeping the point source subtraction because of the paucity of objects in Fatigoni et al. (2021) outside the SRT M31 field. Results are shown in Table 7 (where we show the case from Harper et al. (2023) with no subtraction of sources, as it is the most similar one to our analysis) and the corresponding SEDs in Fig. 9. We recover an AME estimate with reduced significance ($2.2\sigma$) only: $0.93 \pm 0.43$ Jy. Note, though, that the differences between our parameters and those of Harper et al. (2023) are always below the $1\sigma$ level.

We find that the result in our study for the synchrotron index, $\alpha_{syn}$, is consistent with those of Berkhuijsen, Beck & Hoernes (2003) and Planck Collaboration et al. (2015b) ($-0.90 \pm 0.20$ versus $-1.0 \pm 0.2$ and $-0.92 \pm 0.16$, respectively), but steeper than that of Harper et al. (2023) ($-0.63 \pm 0.05$). The synchrotron amplitude, $A_{1 GHz}$ is consistent between this study and Planck Collaboration et al. (2015b) and Harper et al. (2023): $8.6 \pm 1.7$ Jy versus $9.5 \pm 1.1$ and $9.7 \pm 0.6$, respectively. Our free-free estimate is larger, although consistent with a non-detection to within $2\sigma$ too: $5.7 \pm 3.4$ pc cm$^{-6}$ versus $1.8 \pm 1.3$ and $1.1 \pm 0.8$ pc cm$^{-6}$. We get a similar consistency in the AME amplitude with Planck Collaboration et al. (2015b)
<table><caption>**Table 6.** Same as Table 4, but with no point source subtraction applied.</caption><thead><tr><th>Parameter</th><th>Case 1 (all components)</th><th>Case 2 (no AME)</th><th>Case 3 (no free-free)</th><th>From Battistelli et al. (2019)</th></tr></thead><tbody><tr><td>S<sub>1</sub> GHz (Jy)</td><td>7.21 ± 1.15</td><td>7.62 ± 1.13</td><td>8.17 ± 0.98</td><td>6.97 ± 0.54</td></tr><tr><td>α<sub>syn</sub></td><td>-0.88 ± 0.18</td><td>-0.70 ± 0.13</td><td>-0.76 ± 0.16</td><td>-1.1 ± 0.1</td></tr><tr><td>EM (pc cm<sup>-6</sup>)</td><td>4.97 ± 3.11</td><td>7.34 ± 3.84</td><td>—</td><td>—</td></tr><tr><td>A<sub>AME</sub> (Jy)</td><td>0.83 ± 0.37</td><td>—</td><td>1.07 ± 0.44</td><td>—</td></tr><tr><td>ν<sub>AME</sub> (GHz)</td><td>19.49 ± 4.26</td><td>—</td><td>20.58 ± 4.51</td><td>—</td></tr><tr><td>W<sub>AME</sub></td><td>0.58 ± 0.26</td><td>—</td><td>0.74 ± 0.29</td><td>—</td></tr><tr><td>τ<sub>353</sub></td><td>3.76 ± 0.33</td><td>3.79 ± 0.33</td><td>3.73 ± 0.32</td><td>—</td></tr><tr><td>β<sub>d</sub></td><td>1.63 ± 0.09</td><td>1.68 ± 0.08</td><td>1.59 ± 0.08</td><td>1.49 ± 0.056</td></tr><tr><td>T<sub>d</sub> [K]</td><td>19.01 ± 0.84</td><td>18.65 ± 0.78</td><td>19.26 ± 0.83</td><td>18.8 ± 0.54</td></tr><tr><td>S<sup>ff</sup><sub>1 GHz</sub> (Jy)</td><td>0.72 ± 0.45</td><td>1.06 ± 0.56</td><td>—</td><td>0.33 ± 0.26</td></tr><tr><td>S<sup>AME</sup> @ 25 GHz (Jy)</td><td>0.75 ± 0.36</td><td>—</td><td>1.03 ± 0.43</td><td>1.45 ± 0.18</td></tr><tr><td>S<sup>dust</sup><sub>3000 GHz</sub> (Jy)</td><td>3319 ± 1165</td><td>3264 ± 1108</td><td>3356 ± 1129</td><td>3180 ± 230</td></tr><tr><td>χ<sup>red</sup></td><td>0.46</td><td>0.95</td><td>0.36</td><td>0.36</td></tr></tbody></table>

<table><caption>**Table 7.** Same as Table 4, but using an elliptical aperture far from M31 as the background region. Point source subtraction was not applied in this case. This is the scenario most similar to the analysis of Planck Collaboration et al. (2015b). We show the results of Harper et al. (2023) without applying the source subtraction, as it is the most similar case to ours.</caption><thead><tr><th>Parameter</th><th>Case 1 (all components)</th><th>Case 2 (no AME)</th><th>Case 3 (no free-free)</th><th>Planck Collaboration et al. (2015b)</th><th>Harper et al. (2023)</th></tr></thead><tbody><tr><td>S<sub>1</sub> GHz (Jy)</td><td>8.62 ± 1.66</td><td>8.85 ± 1.64</td><td>9.97 ± 1.42</td><td>9.5 ± 1.1</td><td>9.7 ± 0.6</td></tr><tr><td>α<sub>syn</sub></td><td>-0.90 ± 0.20</td><td>-0.75 ± 0.17</td><td>-0.76 ± 0.16</td><td>-0.92 ± 0.16</td><td>-0.63 ± 0.05</td></tr><tr><td>EM (pc cm<sup>-6</sup>)</td><td>5.53 ± 3.29</td><td>8.64 ± 3.88</td><td>—</td><td>1.8 ± 1.3</td><td>1.1 ± 0.8</td></tr><tr><td>A<sub>AME</sub> (Jy)</td><td>0.92 ± 0.43</td><td>—</td><td>1.22 ± 0.52</td><td>—</td><td>—</td></tr><tr><td>ν<sub>AME</sub> (GHz)</td><td>20.84 ± 5.06</td><td>—</td><td>22.46 ± 5.03</td><td>—</td><td>25 ± 2</td></tr><tr><td>τ<sub>353</sub></td><td>3.50 ± 0.32</td><td>3.53 ± 0.31</td><td>3.46 ± 0.30</td><td>—</td><td>—</td></tr><tr><td>β<sub>d</sub></td><td>1.53 ± 0.08</td><td>1.58 ± 0.08</td><td>1.50 ± 0.08</td><td>1.62 ± 0.11</td><td>1.76 ± 0.05</td></tr><tr><td>T<sub>d</sub> [K]</td><td>19.59 ± 0.88</td><td>19.27 ± 0.82</td><td>19.86 ± 0.87</td><td>18.2 ± 1.0</td><td>17.3 ± 0.4</td></tr><tr><td>S<sup>ff</sup><sub>1 GHz</sub> (Jy)</td><td>1.00 ± 0.59</td><td>1.56 ± 0.70</td><td>—</td><td>—</td><td>—</td></tr><tr><td>S<sup>AME</sup> @ 25 GHz (Jy)</td><td>0.74 ± 0.43</td><td>—</td><td>1.13 ± 0.52</td><td>0.7 ± 0.3</td><td>0.32 ± 0.11</td></tr><tr><td>S<sup>dust</sup><sub>3000 GHz</sub> (Jy)</td><td>3933 ± 1358</td><td>3874 ± 1295</td><td>3977 ± 1322</td><td>3340 ± 440</td><td>—</td></tr><tr><td>χ<sup>red</sup></td><td>0.77</td><td>1.12</td><td>0.71</td><td>1.1</td><td>2.6</td></tr></tbody></table>

($S_{\text{AME}}^{30\text{GHz}} = 0.74 \pm 0.43$ Jy versus $S_{\text{AME}}^{30\text{GHz}} = 0.7 \pm 0.3$ Jy), and the value from Harper et al. (2023) ($0.32 \pm 0.11$ Jy) is within $1\sigma$. The AME peak frequency value increased slightly from Section 4 to $\nu_{\text{AME}} = 20.8 \pm 5.1$ GHz, while the AME width decreased to $W_{\text{AME}} = 0.55 \pm 0.29$. Harper et al. (2023) found a higher value for $\nu_{\text{AME}}$, $25 \pm 2$ GHz. However, one must take into account that our model includes an additional free varying parameter driving the width of the AME distribution. On the other hand, both Planck Collaboration et al. (2015b) and Harper et al. (2023) fixed this value by assuming that the AME was fully explained by a single template model (WIM and WNM templates from SPDUST -Ali-Haïmoud, Hirata & Dickinson 2009; Silsbee, Ali-Haïmoud & Hirata 2011-, respectively). Our approach would be more similar to a galaxy hosting multiple AME regions with different spectral dependencies.

Finally, the AME emissivity in this case can be also computed; we find $T_{\text{AME}}^{30\text{GHz}}/\tau_{353} = 4.0 \pm 2.7$ K/$\tau_{353}$, much higher than the one obtained by Harper et al. (2023). The $\epsilon_{\text{AME}}^{28.4\text{GHz}}$ estimate is lower than that from Section 4.2, being $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 8.4 \pm 5.0$ µK MJy$^{-1}$ sr, although within $1\sigma$, because of the large uncertainties.

When evaluating our fit at C-BASS nominal frequency, 4.76 GHz, we find $S_{4.76\text{GHz}} = 3.19 \pm 0.64$ Jy, the large dispersion being consistent with the lack of data between 2–10 GHz. This difference with the value quoted in Harper et al. (2023) ($4.06 \pm 0.14$ Jy) is between the differences for 1.42 GHz and QUIJOTE-MFI surveys (1.5 and 0.55 Jy, respectively). It is noticeable in Fig. 9 that the C-BASS

measurement falls well within the uncertainty of our fitted model. Finally, it is worth noting that C-BASS and QUIJOTE are calibrated to different calibration scales. As explained in Harper et al. (2023), C-BASS data uses the model for Tau A of Weiland et al. (2011). On the contrary, as explained in Rubiño-Martín et al. (2023), QUIJOTE data are calibrated to an updated model that is based on the same data used by Weiland et al. (2011), but including some improvements and the addition of Tau A flux densities from Planck data. Evaluation of the C-BASS model at 4.76 GHz for epoch 2014.4 (C-BASS reference epoch, as observations took place between July 2013 and March 2015) yields a value that is 4.2 per cent higher than the prediction from the QUIJOTE model. Therefore, rescaling the C-BASS point to the QUIJOTE model results in 3.89 Jy, an even better consistency with the fitted models shown in Fig. 9.

## 6 CONCLUSIONS

In this study, we used new data from the QUIJOTE-MFI experiment at 10–20 GHz to study the Andromeda Galaxy (M31). We built the integrated intensity SED of M31 between 0.408 and 3000 GHz and performed a component separation analysis after subtracting the emission from CMB anisotropies and point sources in the region present in the Fatigoni et al. (2021) catalogue. We measure the AME amplitude with a significance of $3.2\sigma$ in the integrated SED of M31, independently of the previous B19 detection. Moreover, our
<table><caption>**Table 6.** Same as Table 4, but with no point source subtraction applied.</caption><thead><tr><th>Parameter</th><th>Case 1 (all components)</th><th>Case 2 (no AME)</th><th>Case 3 (no free-free)</th><th>From Battistelli et al. (2019)</th></tr></thead><tbody><tr><td>S<sub>1</sub> GHz (Jy)</td><td>7.21 ± 1.15</td><td>7.62 ± 1.13</td><td>8.17 ± 0.98</td><td>6.97 ± 0.54</td></tr><tr><td>α<sub>syn</sub></td><td>-0.88 ± 0.18</td><td>-0.70 ± 0.13</td><td>-0.76 ± 0.16</td><td>-1.1 ± 0.1</td></tr><tr><td>EM (pc cm<sup>-6</sup>)</td><td>4.97 ± 3.11</td><td>7.34 ± 3.84</td><td>—</td><td>—</td></tr><tr><td>A<sub>AME</sub> (Jy)</td><td>0.83 ± 0.37</td><td>—</td><td>1.07 ± 0.44</td><td>—</td></tr><tr><td>ν<sub>AME</sub> (GHz)</td><td>19.49 ± 4.26</td><td>—</td><td>20.58 ± 4.51</td><td>—</td></tr><tr><td>W<sub>AME</sub></td><td>0.58 ± 0.26</td><td>—</td><td>0.74 ± 0.29</td><td>—</td></tr><tr><td>τ<sub>353</sub></td><td>3.76 ± 0.33</td><td>3.79 ± 0.33</td><td>3.73 ± 0.32</td><td>—</td></tr><tr><td>β<sub>d</sub></td><td>1.63 ± 0.09</td><td>1.68 ± 0.08</td><td>1.59 ± 0.08</td><td>1.49 ± 0.056</td></tr><tr><td>T<sub>d</sub> [K]</td><td>19.01 ± 0.84</td><td>18.65 ± 0.78</td><td>19.26 ± 0.83</td><td>18.8 ± 0.54</td></tr><tr><td>S<sup>ff</sup><sub>1 GHz</sub> (Jy)</td><td>0.72 ± 0.45</td><td>1.06 ± 0.56</td><td>—</td><td>0.33 ± 0.26</td></tr><tr><td>S<sup>AME</sup> @ 25 GHz (Jy)</td><td>0.75 ± 0.36</td><td>—</td><td>1.03 ± 0.43</td><td>1.45 ± 0.18</td></tr><tr><td>S<sup>dust</sup><sub>3000 GHz</sub> (Jy)</td><td>3319 ± 1165</td><td>3264 ± 1108</td><td>3356 ± 1129</td><td>3180 ± 230</td></tr><tr><td>χ<sup>red</sup></td><td>0.46</td><td>0.95</td><td>0.36</td><td>0.36</td></tr></tbody></table>

<table><caption>**Table 7.** Same as Table 4, but using an elliptical aperture far from M31 as the background region. Point source subtraction was not applied in this case. This is the scenario most similar to the analysis of Planck Collaboration et al. (2015b). We show the results of Harper et al. (2023) without applying the source subtraction, as it is the most similar case to ours.</caption><thead><tr><th>Parameter</th><th>Case 1 (all components)</th><th>Case 2 (no AME)</th><th>Case 3 (no free-free)</th><th>Planck Collaboration et al. (2015b)</th><th>Harper et al. (2023)</th></tr></thead><tbody><tr><td>S<sub>1</sub> GHz (Jy)</td><td>8.62 ± 1.66</td><td>8.85 ± 1.64</td><td>9.97 ± 1.42</td><td>9.5 ± 1.1</td><td>9.7 ± 0.6</td></tr><tr><td>α<sub>syn</sub></td><td>-0.90 ± 0.20</td><td>-0.75 ± 0.17</td><td>-0.76 ± 0.16</td><td>-0.92 ± 0.16</td><td>-0.63 ± 0.05</td></tr><tr><td>EM (pc cm<sup>-6</sup>)</td><td>5.53 ± 3.29</td><td>8.64 ± 3.88</td><td>—</td><td>1.8 ± 1.3</td><td>1.1 ± 0.8</td></tr><tr><td>A<sub>AME</sub> (Jy)</td><td>0.92 ± 0.43</td><td>—</td><td>1.22 ± 0.52</td><td>—</td><td>—</td></tr><tr><td>ν<sub>AME</sub> (GHz)</td><td>20.84 ± 5.06</td><td>—</td><td>22.46 ± 5.03</td><td>—</td><td>25 ± 2</td></tr><tr><td>τ<sub>353</sub></td><td>3.50 ± 0.32</td><td>3.53 ± 0.31</td><td>3.46 ± 0.30</td><td>—</td><td>—</td></tr><tr><td>β<sub>d</sub></td><td>1.53 ± 0.08</td><td>1.58 ± 0.08</td><td>1.50 ± 0.08</td><td>1.62 ± 0.11</td><td>1.76 ± 0.05</td></tr><tr><td>T<sub>d</sub> [K]</td><td>19.59 ± 0.88</td><td>19.27 ± 0.82</td><td>19.86 ± 0.87</td><td>18.2 ± 1.0</td><td>17.3 ± 0.4</td></tr><tr><td>S<sup>ff</sup><sub>1 GHz</sub> (Jy)</td><td>1.00 ± 0.59</td><td>1.56 ± 0.70</td><td>—</td><td>—</td><td>—</td></tr><tr><td>S<sup>AME</sup> @ 25 GHz (Jy)</td><td>0.74 ± 0.43</td><td>—</td><td>1.13 ± 0.52</td><td>0.7 ± 0.3</td><td>0.32 ± 0.11</td></tr><tr><td>S<sup>dust</sup><sub>3000 GHz</sub> (Jy)</td><td>3933 ± 1358</td><td>3874 ± 1295</td><td>3977 ± 1322</td><td>3340 ± 440</td><td>—</td></tr><tr><td>χ<sup>red</sup></td><td>0.77</td><td>1.12</td><td>0.71</td><td>1.1</td><td>2.6</td></tr></tbody></table>

($S_{\text{AME}}^{30\text{GHz}} = 0.74 \pm 0.43$ Jy versus $S_{\text{AME}}^{30\text{GHz}} = 0.7 \pm 0.3$ Jy), and the value from Harper et al. (2023) ($0.32 \pm 0.11$ Jy) is within $1\sigma$. The AME peak frequency value increased slightly from Section 4 to $\nu_{\text{AME}} = 20.8 \pm 5.1$ GHz, while the AME width decreased to $W_{\text{AME}} = 0.55 \pm 0.29$. Harper et al. (2023) found a higher value for $\nu_{\text{AME}}$, $25 \pm 2$ GHz. However, one must take into account that our model includes an additional free varying parameter driving the width of the AME distribution. On the other hand, both Planck Collaboration et al. (2015b) and Harper et al. (2023) fixed this value by assuming that the AME was fully explained by a single template model (WIM and WNM templates from SPDUST -Ali-Haïmoud, Hirata & Dickinson 2009; Silsbee, Ali-Haïmoud & Hirata 2011-, respectively). Our approach would be more similar to a galaxy hosting multiple AME regions with different spectral dependencies.

Finally, the AME emissivity in this case can be also computed; we find $T_{\text{AME}}^{30\text{GHz}}/\tau_{353} = 4.0 \pm 2.7$ K/$(\tau_{353})$, much higher than the one obtained by Harper et al. (2023). The $\epsilon_{\text{AME}}^{28.4\text{GHz}}$ estimate is lower than that from Section 4.2, being $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 8.4 \pm 5.0$ µK MJy$^{-1}$ sr, although within $1\sigma$, because of the large uncertainties.

When evaluating our fit at C-BASS nominal frequency, 4.76 GHz, we find $S_{4.76\text{GHz}} = 3.19 \pm 0.64$ Jy, the large dispersion being consistent with the lack of data between 2–10 GHz. This difference with the value quoted in Harper et al. (2023) ($4.06 \pm 0.14$ Jy) is between the differences for 1.42 GHz and QUIJOTE-MFI surveys (1.5 and 0.55 Jy, respectively). It is noticeable in Fig. 9 that the C-BASS

measurement falls well within the uncertainty of our fitted model. Finally, it is worth noting that C-BASS and QUIJOTE are calibrated to different calibration scales. As explained in Harper et al. (2023), C-BASS data uses the model for Tau A of Weiland et al. (2011). On the contrary, as explained in Rubiño-Martín et al. (2023), QUIJOTE data are calibrated to an updated model that is based on the same data used by Weiland et al. (2011), but including some improvements and the addition of Tau A flux densities from Planck data. Evaluation of the C-BASS model at 4.76 GHz for epoch 2014.4 (C-BASS reference epoch, as observations took place between July 2013 and March 2015) yields a value that is 4.2 per cent higher than the prediction from the QUIJOTE model. Therefore, rescaling the C-BASS point to the QUIJOTE model results in 3.89 Jy, an even better consistency with the fitted models shown in Fig. 9.

## 6 CONCLUSIONS

In this study, we used new data from the QUIJOTE-MFI experiment at 10–20 GHz to study the Andromeda Galaxy (M31). We built the integrated intensity SED of M31 between 0.408 and 3000 GHz and performed a component separation analysis after subtracting the emission from CMB anisotropies and point sources in the region present in the Fatigoni et al. (2021) catalogue. We measure the AME amplitude with a significance of $3.2\sigma$ in the integrated SED of M31, independently of the previous B19 detection. Moreover, our
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Same as Figs 4 and 6, but now using an elliptical aperture far from M31 as the background region. Point source emission has not been subtracted for this case. This case is the one most similar to the analyses of Planck Collaboration et al. (2015b) and Harper et al. (2023). We are adding the C-BASS point from the latter just for comparison purposes.</figcaption></figure>

fitted model is consistent with the 6.6 GHz measurement provided by that study. The addition of QUIJOTE-MFI data improves the definition of the AME spectral parameters: its peak frequency, $v_{AME} = 17.2 \pm 3.2$ GHz, and width, $W_{AME} = 0.58 \pm 0.16$. $v_{AME}$ is low for most theoretical models but consistent with most recent studies pointing to $v_{AME}$ being below 21 GHz (e.g. Harper et al. 2022; Fernández-Torreiro et al. 2023). However, $W_{AME}$ is slightly greater than expected from theoretical models, but again consistent with recent studies (e.g. Fernández-Torreiro et al. 2023; Poidevin et al. 2023). This can be explained by the probable large number of AME regions within M31 with different $v_{AME}$ and $W_{AME}$ values, which

naturally broaden the distribution of the integrated measurement. However, the significance of the AME detection is not as dependent on the QUIJOTE-MFI data addition as on the subtraction of point sources. The significance of the measurement of $A_{AME}$ is still greater than $3\sigma$ (3.1$\sigma$) when discarding the QUIJOTE-MFI data, but this significance falls below $3\sigma$ (2.3$\sigma$) when no point source subtraction is applied. The $A_{AME}$ significance also falls below $3\sigma$ when the COMMANDER map is used to subtract the CMB anisotropies instead of NILC or SMICA, highlighting the importance of the choice of the CMB map. When studying the AME emissivity in M31, we found $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 9.6 \pm 3.1 \text{ µK MJy}^{-1}$ sr compared to the $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 8.8 \pm 3.8 \text{ µK MJy}^{-1}$ sr estimate obtained by Fernández-Torreiro et al. (2023) for the Milky Way, pointing to a similar AME behaviour in both.

We also provide a similar comparison to Harper et al. (2023), who recently measured AME in M31 again, although with much lower amplitude than the previous studies of B19 and of Planck Collaboration et al. (2015b). This difference is driven mainly by the lower flux density derived from SRT data in B19 as compared with the C-BASS value of Harper et al. (2023). However, the two works also show important differences in their data analysis processes. Especially, the different apertures used to perform aperture photometry are of extreme importance. We find consistent results with both when replicating each of the methodologies, while being independent from them, i.e. not using their data. The AME emissivity estimate from Harper et al. (2023) is off by a factor larger than 20 when compared with previous measurements (Planck Collaboration et al. 2015a; Harper et al. 2022; Fernández-Torreiro et al. 2023). When replicating the conditions from Harper et al. (2023), we obtain an AME emissivity estimate $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 8.4 \pm 5.0 \text{ µK MJy}^{-1}$ sr, lower than the previous value of $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 9.6 \pm 3.1 \text{ µK MJy}^{-1}$ sr, probably due to the larger importance of the outskirts of M31. This probably explains why this value is also lower than the estimate for the MW; when observing the MW from its inside, we are not sensitive to these outer regions. This emissivity estimate is equal to $4.0 \pm 2.7 \text{ K/}\tau_{353}$ in $T_{\text{AME}}^{30\text{GHz}}$ units, more in line with previous measurements (Davies et al. 2006; Planck Collaboration et al. 2015a, 2016c; Hensley & Draine 2017; Harper et al. 2022; Fernández-Torreiro et al. 2023) than the one from Harper et al. (2023).

M31 remains the only detection to date of AME in the integrated spectrum of an external galaxy. We have also provided upper limits for the polarization fraction of this AME component: 30 and 20 per cent at 11 and 13 GHz, respectively. Although these limits are large when compared with those from brighter AME sources within our Galaxy, this is the first time that such upper limits have been computed for an extragalactic object. The results of this work encourage the reproduction of this kind of measurements for galaxies with different sets of properties throughout the lifetime of the Universe.

## ACKNOWLEDGEMENTS

We thank the anonymous referee whose comments helped to improve this work. We thank the staff of the Teide Observatory for invaluable assistance in the commissioning and operation of QUIJOTE. The QUIJOTE experiment is being developed by the Instituto de Astrofísica de Canarias (IAC), the Instituto de Física de Cantabria (IFCA), and the Universities of Cantabria, Manchester and Cambridge. Partial financial support was provided by the Spanish Ministry of Science and Innovation under the projects AYA2007-68058-C03-01, AYA2007-
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Same as Figs 4 and 6, but now using an elliptical aperture far from M31 as the background region. Point source emission has not been subtracted for this case. This case is the one most similar to the analyses of Planck Collaboration et al. (2015b) and Harper et al. (2023). We are adding the C-BASS point from the latter just for comparison purposes.</figcaption></figure>

fitted model is consistent with the 6.6 GHz measurement provided by that study. The addition of QUIJOTE-MFI data improves the definition of the AME spectral parameters: its peak frequency, $v_{AME} = 17.2 \pm 3.2$ GHz, and width, $W_{AME} = 0.58 \pm 0.16$. $v_{AME}$ is low for most theoretical models but consistent with most recent studies pointing to $v_{AME}$ being below 21 GHz (e.g. Harper et al. 2022; Fernández-Torreiro et al. 2023). However, $W_{AME}$ is slightly greater than expected from theoretical models, but again consistent with recent studies (e.g. Fernández-Torreiro et al. 2023; Poidevin et al. 2023). This can be explained by the probable large number of AME regions within M31 with different $v_{AME}$ and $W_{AME}$ values, which

naturally broaden the distribution of the integrated measurement. However, the significance of the AME detection is not as dependent on the QUIJOTE-MFI data addition as on the subtraction of point sources. The significance of the measurement of $A_{AME}$ is still greater than $3\sigma$ (3.1$\sigma$) when discarding the QUIJOTE-MFI data, but this significance falls below $3\sigma$ (2.3$\sigma$) when no point source subtraction is applied. The $A_{AME}$ significance also falls below $3\sigma$ when the COMMANDER map is used to subtract the CMB anisotropies instead of NILC or SMICA, highlighting the importance of the choice of the CMB map. When studying the AME emissivity in M31, we found $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 9.6 \pm 3.1 \text{ µK MJy}^{-1}$ sr compared to the $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 8.8 \pm 3.8 \text{ µK MJy}^{-1}$ sr estimate obtained by Fernández-Torreiro et al. (2023) for the Milky Way, pointing to a similar AME behaviour in both.

We also provide a similar comparison to Harper et al. (2023), who recently measured AME in M31 again, although with much lower amplitude than the previous studies of B19 and of Planck Collaboration et al. (2015b). This difference is driven mainly by the lower flux density derived from SRT data in B19 as compared with the C-BASS value of Harper et al. (2023). However, the two works also show important differences in their data analysis processes. Especially, the different apertures used to perform aperture photometry are of extreme importance. We find consistent results with both when replicating each of the methodologies, while being independent from them, i.e. not using their data. The AME emissivity estimate from Harper et al. (2023) is off by a factor larger than 20 when compared with previous measurements (Planck Collaboration et al. 2015a; Harper et al. 2022; Fernández-Torreiro et al. 2023). When replicating the conditions from Harper et al. (2023), we obtain an AME emissivity estimate $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 8.4 \pm 5.0 \text{ µK MJy}^{-1}$ sr, lower than the previous value of $\epsilon_{\text{AME}}^{28.4\text{GHz}} = 9.6 \pm 3.1 \text{ µK MJy}^{-1}$ sr, probably due to the larger importance of the outskirts of M31. This probably explains why this value is also lower than the estimate for the MW; when observing the MW from its inside, we are not sensitive to these outer regions. This emissivity estimate is equal to $4.0 \pm 2.7 \text{ K/}\tau_{353}$ in $T_{\text{AME}}^{30\text{GHz}}$ units, more in line with previous measurements (Davies et al. 2006; Planck Collaboration et al. 2015a, 2016c; Hensley & Draine 2017; Harper et al. 2022; Fernández-Torreiro et al. 2023) than the one from Harper et al. (2023).

M31 remains the only detection to date of AME in the integrated spectrum of an external galaxy. We have also provided upper limits for the polarization fraction of this AME component: 30 and 20 per cent at 11 and 13 GHz, respectively. Although these limits are large when compared with those from brighter AME sources within our Galaxy, this is the first time that such upper limits have been computed for an extragalactic object. The results of this work encourage the reproduction of this kind of measurements for galaxies with different sets of properties throughout the lifetime of the Universe.

## ACKNOWLEDGEMENTS

We thank the anonymous referee whose comments helped to improve this work. We thank the staff of the Teide Observatory for invaluable assistance in the commissioning and operation of QUIJOTE. The QUIJOTE experiment is being developed by the Instituto de Astrofísica de Canarias (IAC), the Instituto de Física de Cantabria (IFCA), and the Universities of Cantabria, Manchester and Cambridge. Partial financial support was provided by the Spanish Ministry of Science and Innovation under the projects AYA2007-68058-C03-01, AYA2007-
68058-C03-02, AYA2010-21766-C03-01, AYA2010-21766-C03-02, AYA2014-60438-P, ESP2015-70646-C2-1-R, AYA2017-84185-P, ESP2017-83921-C2-1-R, PID2019-110610RB-C21, PID2020-120514GB-I00, IACA13-3E-2336, IACA15-BE-3707, EQC2018-004918-P, the Severo Ochoa Programs SEV-2015-0548 and CEX2019-000920-S, the Maria de Maeztu Program MDM-2017-0765, and by the Consolider-Ingenio project CSD2010-00064 (EPI: Exploring the Physics of Inflation). We acknowledge support from the ACIISI, Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias, and the European Regional Development Fund (ERDF) under grant with reference ProID2020010108. This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement number 687312 (RADIOFOREGROUNDS). MFT acknowledges support from the Spanish Agencia Estatal de Investigación (AEI) of the Ministerio de Ciencia, Innovación y Universidades (MCIU) and the European Social Fund (ESF) under grant with reference PRE-C-2018-0067. CA-T acknowledges support from the Millennium Nucleus on Young Exoplanets and their Moons (YEMS). FP acknowledges support from the Agencia Canaria de Investigación, Innovación y Sociedad de la Información (ACIISI) under the European FEDER (Fondo Europeo de Desarrollo Regional) de Canarias 2014–2020 grant No. PROID2021010078.

We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center. We acknowledge the use of data provided by the Centre d'Analyse de Données Etendues (CADE), a service of IRAP-UPS/CNRS (http://cade.irap.omp.eu, Paradis et al. 2012). This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France (Wenger et al. 2000). Based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. Some of the presented results are based on observations obtained with the QUIJOTE experiment (http://rese.arch.iac.es/proyecto/quijote). Some of the results in this paper have been derived using the healpy and HEALPix packages (Górski et al. 2005; Zonca et al. 2019). We have also used SCIPY (Virtanen et al. 2020), EMCEE (Foreman-Mackey et al. 2013), NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), CORNER (Foreman-Mackey 2016) and ASTROPY (Astropy Collaboration et al. 2013, 2018) PYTHON packages.

## DATA AVAILABILITY

The QUIJOTE raster scan data used in this paper are property of the QUIJOTE Collaboration and can only be shared on request to the corresponding authors. All M31 fits are available upon request to the QUIJOTE collaboration. QUIJOTE-MFI Wide Survey data is available in the QUIJOTE collaboration webpage (https://rese.arch.iac.es/proyecto/quijote).

## REFERENCES

* Akaike H., 1974, IEEE Trans. Autom. Control, 19, 716

* Ali-Haïmoud Y., Hirata C. M., Dickinson C., 2009, MNRAS, 395, 1055

* Astropy Collaboration et al., 2013, A&A, 558, A33

* Astropy Collaboration et al., 2018, AJ, 156, 123

* Battistelli E. S. et al., 2019, ApJ, 877, L31 (B19)

* Bennett C. L. et al., 2013, ApJS, 208, 20

* Berkhuijsen E. M., 1972, A&AS, 5, 263

* Berkhuijsen E. M., Beck R., Hoernes P., 2003, A&A, 398, 937

* Bianchi S. et al., 2022, A&A, 658, L8

* Cepeda-Arroita R. et al., 2021, MNRAS, 503, 2927

* Davies R. D., Dickinson C., Banday A. J., Jaffe T. R., Górski K. M., Davis R. J., 2006, MNRAS, 370, 1125

* de Vaucouleurs G., de Vaucouleurs A., Corwin Herold G. J., Buta R. J., Paturel G., Fouque P., 1991, Third Reference Catalogue of Bright Galaxies. Springer, New York, NY (USA)

* Dickinson C. et al., 2018, New A Rev., 80, 1

* Dickinson C., Davies R. D., Davis R. J., 2003, MNRAS, 341, 369

* Dickinson C., Peel M., Vidal M., 2011, MNRAS, 418, L35

* Draine B. T., Lazarian A., 1998a, ApJ, 494, L19

* Draine B. T., Lazarian A., 1998b, ApJ, 508, 157

* Fatigoni S. et al., 2021, A&A, 651, A98

* Fernández-Torreiro M. et al., 2023, MNRAS 526 1343

* Finkbeiner D. P., 2003, ApJS, 146, 407

* Ford G. P. et al., 2013, ApJ, 769, 55

* Foreman-Mackey D., 2016, J. Open Source Softw., 1, 24

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Fritz J. et al., 2012, A&A, 546, A34

* Génova-Santos R. et al., 2015, MNRAS, 452, 4169

* Génova-Santos R. et al., 2017, MNRAS, 464, 4107

* Goodwin S. P., Gribbin J., Hendry M. A., 1998, The Observatory, 118, 201

* Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

* Guidi F. et al., 2021, MNRAS, 507, 3707

* Guidi F. et al., 2023, MNRAS, 519, 3460

* Harper S. E. et al., 2022, MNRAS, 513, 5900

* Harper S. E. et al., 2023, MNRAS, 523, 3471

* Harris C. R. et al., 2020, Nature, 585, 357

* Haslam C. G. T., Salter C. J., Stoffel H., Wilson W. E., 1982, A&AS, 47, 1

* Hauser M. G. et al., 1998, ApJ, 508, 25

* Healey S. E., Romani R. W., Taylor G. B., Sadler E. M., Ricci R., Murphy T., Ulvestad J. S., Winn J. N., 2007, ApJS, 171, 61

* Hensley B. S., Draine B. T., 2017, ApJ, 836, 179

* Hensley B., Murphy E., Staguhn J., 2015, MNRAS, 449, 809

* Herman D. et al., 2023, A&A, 675, A15

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

* Kogut A., Banday A. J., Bennett C. L., Gorski K. M., Hinshaw G., Smoot G. F., Wright E. I., 1996, ApJ, 464, L5

* Leitch E. M., Readhead A. C. S., Pearson T. J., Myers S. T., 1997, ApJ, 486, L23

* López-Caraballo C. H. et al., 2023, preprint (arXiv:2307.15518)

* López-Caraballo C. H., Rubiño-Martín J. A., Rebolo R., Génova-Santos R., 2011, ApJ, 729, 25

* Macellari N., Pierpaoli E., Dickinson C., Vaillancourt J. E., 2011, MNRAS, 418, 888

* McConnachie A. W., Irwin M. J., Ferguson A. M. N., Ibata R. A., Lewis G. F., Tanvir N., 2005, MNRAS, 356, 979

* Murphy E. J. et al., 2010, ApJ, 709, L108

* Murphy E. J. et al., 2012, ApJ, 761, 97

* Murphy E. J., Linden S. T., Dong D., Hensley B. S., Momjian E., Helou G., Evans A. S., 2018, ApJ, 862, 20

* Paradis D., Dobashi K., Shimoikura T., Kawamura A., Onishi T., Fukui Y., Bernard J. P., 2012, A&A, 543, A103

* Peel M. W., Dickinson C., Davies R. D., Clements D. L., Beswick R. J., 2011, MNRAS, 416, L99

* Peel M. W., Genova-Santos R., Dickinson C., Leahy J. P., López-Caraballo C., Fernández-Torreiro M., Rubiño-Martín J. A., Spencer L. D., 2022, RNAAS, 6, 252

* Planck Collaboration I, 2020a, A&A, 641, A1

* Planck Collaboration IV, 2020b, A&A, 641, A4

* Planck Collaboration X, 2016b, A&A, 594, A10
68058-C03-02, AYA2010-21766-C03-01, AYA2010-21766-C03-02, AYA2014-60438-P, ESP2015-70646-C2-1-R, AYA2017-84185-P, ESP2017-83921-C2-1-R, PID2019-110610RB-C21, PID2020-120514GB-I00, IACA13-3E-2336, IACA15-BE-3707, EQC2018-004918-P, the Severo Ochoa Programs SEV-2015-0548 and CEX2019-000920-S, the Maria de Maeztu Program MDM-2017-0765, and by the Consolider-Ingenio project CSD2010-00064 (EPI: Exploring the Physics of Inflation). We acknowledge support from the ACIISI, Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias, and the European Regional Development Fund (ERDF) under grant with reference ProID2020010108. This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement number 687312 (RADIOFOREGROUNDS). MFT acknowledges support from the Spanish Agencia Estatal de Investigación (AEI) of the Ministerio de Ciencia, Innovación y Universidades (MCIU) and the European Social Fund (ESF) under grant with reference PRE-C-2018-0067. CA-T acknowledges support from the Millennium Nucleus on Young Exoplanets and their Moons (YEMS). FP acknowledges support from the Agencia Canaria de Investigación, Innovación y Sociedad de la Información (ACIISI) under the European FEDER (Fondo Europeo de Desarrollo Regional) de Canarias 2014–2020 grant No. PROID2021010078.

We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center. We acknowledge the use of data provided by the Centre d'Analyse de Données Etendues (CADE), a service of IRAP-UPS/CNRS (http://cade.irap.omp.eu, Paradis et al. 2012). This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France (Wenger et al. 2000). Based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. Some of the presented results are based on observations obtained with the QUIJOTE experiment (http://rese arch.iac.es/proyecto/quijote). Some of the results in this paper have been derived using the healpy and HEALPix packages (Górski et al. 2005; Zonca et al. 2019). We have also used SCIPY (Virtanen et al. 2020), EMCEE (Foreman-Mackey et al. 2013), NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), CORNER (Foreman-Mackey 2016) and ASTROPY (Astropy Collaboration et al. 2013, 2018) PYTHON packages.

## DATA AVAILABILITY

The QUIJOTE raster scan data used in this paper are property of the QUIJOTE Collaboration and can only be shared on request to the corresponding authors. All M31 fits are available upon request to the QUIJOTE collaboration. QUIJOTE-MFI Wide Survey data is available in the QUIJOTE collaboration webpage (https://rese arch.iac.es/proyecto/quijote).

## REFERENCES

* Akaike H., 1974, IEEE Trans. Autom. Control, 19, 716

* Ali-Haïmoud Y., Hirata C. M., Dickinson C., 2009, MNRAS, 395, 1055

* Astropy Collaboration et al., 2013, A&A, 558, A33

* Astropy Collaboration et al., 2018, AJ, 156, 123

* Battistelli E. S. et al., 2019, ApJ, 877, L31 (B19)

* Bennett C. L. et al., 2013, ApJS, 208, 20

* Berkhuijsen E. M., 1972, A&AS, 5, 263

* Berkhuijsen E. M., Beck R., Hoernes P., 2003, A&A, 398, 937

* Bianchi S. et al., 2022, A&A, 658, L8

* Cepeda-Arroita R. et al., 2021, MNRAS, 503, 2927

* Davies R. D., Dickinson C., Banday A. J., Jaffe T. R., Górski K. M., Davis R. J., 2006, MNRAS, 370, 1125

* de Vaucouleurs G., de Vaucouleurs A., Corwin Herold G. J., Buta R. J., Paturel G., Fouque P., 1991, Third Reference Catalogue of Bright Galaxies. Springer, New York, NY (USA)

* Dickinson C. et al., 2018, New A Rev., 80, 1

* Dickinson C., Davies R. D., Davis R. J., 2003, MNRAS, 341, 369

* Dickinson C., Peel M., Vidal M., 2011, MNRAS, 418, L35

* Draine B. T., Lazarian A., 1998a, ApJ, 494, L19

* Draine B. T., Lazarian A., 1998b, ApJ, 508, 157

* Fatigoni S. et al., 2021, A&A, 651, A98

* Fernández-Torreiro M. et al., 2023, MNRAS 526 1343

* Finkbeiner D. P., 2003, ApJS, 146, 407

* Ford G. P. et al., 2013, ApJ, 769, 55

* Foreman-Mackey D., 2016, J. Open Source Softw., 1, 24

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Fritz J. et al., 2012, A&A, 546, A34

* Génova-Santos R. et al., 2015, MNRAS, 452, 4169

* Génova-Santos R. et al., 2017, MNRAS, 464, 4107

* Goodwin S. P., Gribbin J., Hendry M. A., 1998, The Observatory, 118, 201

* Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

* Guidi F. et al., 2021, MNRAS, 507, 3707

* Guidi F. et al., 2023, MNRAS, 519, 3460

* Harper S. E. et al., 2022, MNRAS, 513, 5900

* Harper S. E. et al., 2023, MNRAS, 523, 3471

* Harris C. R. et al., 2020, Nature, 585, 357

* Haslam C. G. T., Salter C. J., Stoffel H., Wilson W. E., 1982, A&AS, 47, 1

* Hauser M. G. et al., 1998, ApJ, 508, 25

* Healey S. E., Romani R. W., Taylor G. B., Sadler E. M., Ricci R., Murphy T., Ulvestad J. S., Winn J. N., 2007, ApJS, 171, 61

* Hensley B. S., Draine B. T., 2017, ApJ, 836, 179

* Hensley B., Murphy E., Staguhn J., 2015, MNRAS, 449, 809

* Herman D. et al., 2023, A&A, 675, A15

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

* Kogut A., Banday A. J., Bennett C. L., Gorski K. M., Hinshaw G., Smoot G. F., Wright E. I., 1996, ApJ, 464, L5

* Leitch E. M., Readhead A. C. S., Pearson T. J., Myers S. T., 1997, ApJ, 486, L23

* López-Caraballo C. H. et al., 2023, preprint (arXiv:2307.15518)

* López-Caraballo C. H., Rubiño-Martín J. A., Rebolo R., Génova-Santos R., 2011, ApJ, 729, 25

* Macellari N., Pierpaoli E., Dickinson C., Vaillancourt J. E., 2011, MNRAS, 418, 888

* McConnachie A. W., Irwin M. J., Ferguson A. M. N., Ibata R. A., Lewis G. F., Tanvir N., 2005, MNRAS, 356, 979

* Murphy E. J. et al., 2010, ApJ, 709, L108

* Murphy E. J. et al., 2012, ApJ, 761, 97

* Murphy E. J., Linden S. T., Dong D., Hensley B. S., Momjian E., Helou G., Evans A. S., 2018, ApJ, 862, 20

* Paradis D., Dobashi K., Shimoikura T., Kawamura A., Onishi T., Fukui Y., Bernard J. P., 2012, A&A, 543, A103

* Peel M. W., Dickinson C., Davies R. D., Clements D. L., Beswick R. J., 2011, MNRAS, 416, L99

* Peel M. W., Genova-Santos R., Dickinson C., Leahy J. P., López-Caraballo C., Fernández-Torreiro M., Rubiño-Martín J. A., Spencer L. D., 2022, RNAAS, 6, 252

* Planck Collaboration I, 2020a, A&A, 641, A1

* Planck Collaboration IV, 2020b, A&A, 641, A4

* Planck Collaboration X, 2016b, A&A, 594, A10
Planck Collaboration XI, 2014b, A&A, 571, A11

Planck Collaboration XV, 2014a, A&A, 565, A103

Planck Collaboration XXIII, 2015a, A&A, 580, A13

Planck Collaboration XXIX, 2016a, A&A, 586, A132

Planck Collaboration XXV, 2015b, A&A, 582, A28

Planck Collaboration XXV, 2016c, A&A, 594, A25

Poidevin F. et al., 2019, MNRAS, 486, 462

Poidevin F. et al., 2023, MNRAS, 519, 3481

Reich P., Reich W., 1988, A&AS, 74, 7

Reich P., Testori J. C., Reich W., 2001, A&A, 376, 861

Reich W., 1982, A&AS, 48, 219

Remazeilles M., Dickinson C., Banday A. J., Bigot-Sazy M. A., Ghosh T., 2015, MNRAS, 451, 4311

Richards J. L. et al., 2011, ApJS, 194, 29

Rubiño-Martín J. A. et al., 2010, Ap&SS Proc., 14, 127

Rubiño-Martín J. A. et al., 2023, MNRAS, 519, 3383

Rubiño-Martín J. A., López-Caraballo C. H., Génova-Santos R., Rebolo R., 2012, Adv. Astron., 2012, 351836

Rys S., Machalski J., 1990, A&A, 236, 15

Schwarz G., 1978, Ann. Stat., 6, 461

Silsbee K., Ali-Haïmoud Y., Hirata C. M., 2011, MNRAS, 411, 2750

Stevenson M. A., 2014, ApJ, 781, 113

Tibbs C. T., Israel F. P., Laureijs R. J., Tauber J. A., Partridge B., Peel M. W., Fauvet L., 2018, MNRAS, 477, 4968

Tibbs C. T., Paladini R., Dickinson C., 2012, Adv. Astron., 2012, 124931

Tramonte D. et al., 2023, MNRAS, 519, 3432

Vaillancourt J. E., 2006, PASP, 118, 1340

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Watson R. A., Rebolo R., Rubiño-Martín J. A., Hildebrandt S., Gutiérrez C. M., Fernández-Cerezo S., Hoyland R. J., Battistelli E. S., 2005, ApJ, 624, L89

Weiland J. L. et al., 2011, ApJS, 192, 19

Wenger M. et al., 2000, A&AS, 143, 9

Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, J. Open Source Softw., 4, 1298

This paper has been typeset from a TeX/LA TeX file prepared by the author.

<footnote n="1">¹Instituto de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="2">²Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="3">³Imperial College London, Blackett Lab, Prince Consort Road, London SW7 2AZ, UK</footnote>

<footnote n="4">⁴Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="5">⁵Millennium Nucleus on Young Exoplanets and their Moons (YEMS), 9170124, Santiago Chile</footnote>

<footnote n="6">⁶Consejo Superior de Investigaciones Científicas, E-28006 Madrid, Spain</footnote>

<footnote n="7">⁷Universidad de Cantabria, Departamento de Ingeniería de Comunicaciones, Edificio Ingenieria de Telecomunicación, Plaza de la Ciencia nº 1, E-39005 Santander, Spain</footnote>

<footnote n="8">⁸Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="9">⁹Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="10">¹⁰Instituto de Física de Cantabria (IFCA), CSIC-Univ. de Cantabria, Avda. los Castros, s/n, E-39005 Santander, Spain</footnote>

<footnote n="11">¹¹Departamento de Física Moderna, Universidad de Cantabria, Avda. de los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote n="12">¹²CNRS-UCB International Research Laboratory, Centre Pierre Binétruy, IRL2007, CPB-IN2P3, Berkeley, CA 94720, USA</footnote>

<footnote n="13">¹³Institut d'Astrophysique de Paris, UMR 7095, CNRS & Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="14">¹⁴Jodrell Bank Centre for Astrophysics, Alan Turing Building, Department of Physics & Astronomy, School of Natural Sciences, The University of Manchester, Oxford Road, Manchester, M13 9PL, U.K.</footnote>

<footnote n="15">¹⁵Departamento de Física. Facultad de Ciencias. Universidad de Córdoba. Campus de Rabanales, Edif. C2. Planta Baja. E-14071 Córdoba, Spain</footnote>

<footnote n="16">¹⁶Department of Physics, Xi'an Jiaotong-Liverpool University, 111 Ren'ai Road, Suzhou Dushu Lake Science and Education Innovation District, Suzhou Industrial Park, Suzhou 215123, P.R. China</footnote>
Planck Collaboration XI, 2014b, A&A, 571, A11

Planck Collaboration XV, 2014a, A&A, 565, A103

Planck Collaboration XXIII, 2015a, A&A, 580, A13

Planck Collaboration XXIX, 2016a, A&A, 586, A132

Planck Collaboration XXV, 2015b, A&A, 582, A28

Planck Collaboration XXV, 2016c, A&A, 594, A25

Poidevin F. et al., 2019, MNRAS, 486, 462

Poidevin F. et al., 2023, MNRAS, 519, 3481

Reich P., Reich W., 1988, A&AS, 74, 7

Reich P., Testori J. C., Reich W., 2001, A&A, 376, 861

Reich W., 1982, A&AS, 48, 219

Remazeilles M., Dickinson C., Banday A. J., Bigot-Sazy M. A., Ghosh T., 2015, MNRAS, 451, 4311

Richards J. L. et al., 2011, ApJS, 194, 29

Rubiño-Martín J. A. et al., 2010, Ap&SS Proc., 14, 127

Rubiño-Martín J. A. et al., 2023, MNRAS, 519, 3383

Rubiño-Martín J. A., López-Caraballo C. H., Génova-Santos R., Rebolo R., 2012, Adv. Astron., 2012, 351836

Rys S., Machalski J., 1990, A&A, 236, 15

Schwarz G., 1978, Ann. Stat., 6, 461

Silsbee K., Ali-Haïmoud Y., Hirata C. M., 2011, MNRAS, 411, 2750

Stevenson M. A., 2014, ApJ, 781, 113

Tibbs C. T., Israel F. P., Laureijs R. J., Tauber J. A., Partridge B., Peel M. W., Fauvet L., 2018, MNRAS, 477, 4968

Tibbs C. T., Paladini R., Dickinson C., 2012, Adv. Astron., 2012, 124931

Tramonte D. et al., 2023, MNRAS, 519, 3432

Vaillancourt J. E., 2006, PASP, 118, 1340

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Watson R. A., Rebolo R., Rubiño-Martín J. A., Hildebrandt S., Gutiérrez C. M., Fernández-Cerezo S., Hoyland R. J., Battistelli E. S., 2005, ApJ, 624, L89

Weiland J. L. et al., 2011, ApJS, 192, 19

Wenger M. et al., 2000, A&AS, 143, 9

Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, J. Open Source Softw., 4, 1298

This paper has been typeset from a TeX/LA TeX file prepared by the author.

<footnote n="1">¹Instituto de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="2">²Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="3">³Imperial College London, Blackett Lab, Prince Consort Road, London SW7 2AZ, UK</footnote>

<footnote n="4">⁴Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="5">⁵Millennium Nucleus on Young Exoplanets and their Moons (YEMS), 9170124, Santiago Chile</footnote>

<footnote n="6">⁶Consejo Superior de Investigaciones Científicas, E-28006 Madrid, Spain</footnote>

<footnote n="7">⁷Universidad de Cantabria, Departamento de Ingeniería de Comunicaciones, Edificio Ingenieria de Telecomunicación, Plaza de la Ciencia nº 1, E-39005 Santander, Spain</footnote>

<footnote n="8">⁸Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="9">⁹Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="10">¹⁰Instituto de Física de Cantabria (IFCA), CSIC-Univ. de Cantabria, Avda. los Castros, s/n, E-39005 Santander, Spain</footnote>

<footnote n="11">¹¹Departamento de Física Moderna, Universidad de Cantabria, Avda. de los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote n="12">¹²CNRS-UCB International Research Laboratory, Centre Pierre Binétruy, IRL2007, CPB-IN2P3, Berkeley, CA 94720, USA</footnote>

<footnote n="13">¹³Institut d'Astrophysique de Paris, UMR 7095, CNRS & Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="14">¹⁴Jodrell Bank Centre for Astrophysics, Alan Turing Building, Department of Physics & Astronomy, School of Natural Sciences, The University of Manchester, Oxford Road, Manchester, M13 9PL, U.K.</footnote>

<footnote n="15">¹⁵Departamento de Física. Facultad de Ciencias. Universidad de Córdoba. Campus de Rabanales, Edif. C2. Planta Baja. E-14071 Córdoba, Spain</footnote>

<footnote n="16">¹⁶Department of Physics, Xi'an Jiaotong-Liverpool University, 111 Ren'ai Road, Suzhou Dushu Lake Science and Education Innovation District, Suzhou Industrial Park, Suzhou 215123, P.R. China</footnote>