The structure of the stellar halo of the Andromeda galaxy explored with the NB515 for Subaru/HSC – I. New insights on the stellar halo up to 120 kpc Open Access

ABSTRACT

We analyse the M31 halo and its substructure within a projected radius of 120 kpc using a combination of Subaru/HSC NB515 and Canada France Hawaii Telescope/MegaCam g and i bands. We succeed in separating M31’s halo stars from foreground contamination with |$\sim$|90 per cent accuracy by using the surface gravity sensitive NB515 filter. Based on the selected M31 halo stars, we discover three new substructures, which associate with the Giant Southern Stream (GSS) based on their photometric metallicity estimates. We also produce the distance and photometric metallicity estimates for the known substructures. While these quantities for the GSS are reproduced in our study, we find that the north-western stream shows a steeper distance gradient than found in an earlier study, suggesting that it is likely to have formed in an orbit closer to the Milky Way. For two streams in the eastern halo (Stream C and D), we identify distance gradients that had not been resolved. Finally, we investigate the global halo photometric metallicity distribution and surface brightness profile using the NB515-selected halo stars. We find that the surface brightness of the metal-poor and metal-rich halo populations, and the all population can be fitted to a power-law profile with an index of |$\alpha =-1.65\pm 0.02$|⁠, |$-2.82\pm 0.01$|⁠, and |$-2.44\pm 0.01$|⁠, respectively. In contrast to the relative smoothness of the halo profile, its photometric metallicity distribution appears to be spatially non-uniform with non-monotonic trends with radius, suggesting that the halo population had insufficient time to dynamically homogenize the accreted populations.

1 INTRODUCTION

Understanding how galaxies have formed and evolved to what we see today is one of the most critical issues in astronomy. Searle & Zinn (1978) proposed that large galaxies such as the Milky Way (MW) and the Andromeda Galaxy (M31) have been formed by the accretions of small stellar systems over a long time, based on the lack of spatial dependence in the metallicity of globular clusters. This galaxy formation scenario is consistent with the Lambda-dominated cold dark matter (⁠|$\Lambda$|CDM) theory. The hierarchical structure formation scenario based on |$\Lambda$|CDM model is also suitable for describing large-scale structures beyond the galactic scale (⁠|$\sim$|1 Mpc) and is therefore widely accepted. This model predicts that the accreted stellar systems are disrupted by the tidal forces of the host galaxy, leaving behind tidal streams (Bullock & Johnston 2005; Johnston et al. 2008; Font et al. 2011).

Stellar haloes are a key to understanding galaxy formation because they retain past chemodynamical information. Many of the halo stars found in the MW have low metallicities and high-velocity dispersions, which differ from the properties of disc stars (e.g. Feltzing & Chiba 2013). Analyses of the kinematics and chemical abundances in Galactic halo stars have uncovered that the Galactic halo consists of at least two distinct components, a flattened inner halo with a higher mean metallicity at |$R_{\rm gal} \lt 10~{\rm kpc}$| and a spherical outer halo with a lower mean metallicity at |$R_{\rm gal} \gt 30~{\rm kpc}$| (Carollo et al. 2007, 2010; Conroy et al. 2019). Thanks to the Gaia satellite mission, it has become evident that a large fraction of the field halo stars in the solar neighbourhood exhibit large orbital eccentricities and lower |$[\mathrm{\alpha /Fe}]$| ratios when compared to the majority of stars with similar metallicities. These results are currently interpreted as a stellar population that was accreted into the MW approximately 10 billion years ago (Belokurov et al. 2018; Helmi et al. 2018; Gallart et al. 2019). In addition to these halo stellar populations, observations with the Sloan Digital Sky Survey (York et al. 2000), the Gaia satellite (Gaia Collaboration 2016), and others have discovered over 70 substructures (e.g. Ibata, Gilmore & Irwin 1994; Belokurov et al. 2006, 2018; Mateu et al. 2017; Mateu, Read & Kawata 2018; Shipp et al. 2018). These substructures are evidence that stellar haloes are formed by the accretion of small stellar systems. Furthermore, numerical simulations of the orbits of accreted dwarf galaxies based on the |$\Lambda$|CDM theory revealed that past accretion events have caused dynamical heating of the Galactic disc (e.g. Villalobos & Helmi 2008; Cooper et al. 2010; Grand et al. 2017) and that in situ star formation was predicted in the Galactic halo (e.g. Zolotov et al. 2009, 2010), such in situ formation was observationally detected (e.g. Gallart et al. 2019; Belokurov & Kravtsov 2022). Thus, a detailed study of Galactic stellar haloes is an important clue to understanding galaxy formation.

Similar to the MW, the resolved disc and halo stars in our neighbour, M31, have also been observed and studied in detail. In particular, its low disc inclination of |$12^{\circ }5$| (Simien et al. 1978) allows M31 to provide a comprehensive view of the halo. Therefore, M31 is an excellent system for studying galaxy formation.

Many photometric and spectroscopic observations in the M31 halo have been conducted, including large-scale surveys such as Pan-Andromeda Archaeological Survey (PAndAS; McConnachie et al. 2009, 2018) with the Canada France Hawaii Telescope (CFHT)/MegaCam and the Spectroscopic and Photometric Landscape of Andromeda’s Stellar Halo (SPLASH; Gilbert et al. 2009) survey with the Keck II/DEep Imaging Multi-Object Spectrograph (DEIMOS). The results of these observations show that there are some similarities to the MW halo, but various differences have also been found. Like the Galactic halo, the surface brightness of the M31 outer halo follows a power-law profile, and many stars have low metallicities (⁠|$[\mathrm{Fe/H}]\lt -1.5$|⁠) (Guhathakurta et al. 2005; Kalirai et al. 2006; Gilbert et al. 2014, 2020; Ibata et al. 2014 ). On the other hand, the inner halo of M31 has been found to have higher metallicities (⁠|$[\mathrm{Fe/H}]\sim -0.6$|⁠) and younger ages (⁠|$\sim 8 \,\mathrm{Gyr}$|⁠) than the Galactic inner halo (e.g. Brown et al. 2003, 2006; Gilbert et al. 2014). In addition to the differences in the metallicities, the masses of stellar haloes of the MW and M31 are very different; the stellar halo of the MW has |$\sim 10^9~{\rm M}_{\odot }$| (Deason, Belokurov & Sanders 2019), and the accreted masses of M31 was expected to be |$\gt 1-2\times 10^{10}~{\rm M}_{\odot }$| (D’Souza & Bell 2018).

Substructures have also been discovered in the M31 halo similar to the Galactic halo. Starting with the discovery of the Giant Southern Stream (GSS) by Ibata et al. (2001), more than 10 substructures have been discovered in the eastern and western parts of the M31 disc and in the north-western and south-western parts of the halo (e.g. Ibata et al. 2001; Ferguson et al. 2002, 2005; Ibata et al. 2007; McConnachie et al. 2009; Richardson et al. 2011).

As in the Galaxy, these observations provide many opportunities for simulations. In particular, GSS formation scenarios have been put forward through numerical simulations: this substructure can be formed from radial merger with a dwarf galaxy of masses of |$\sim 10^9$||$\mathrm{M}_{\odot }~$| (e.g. Fardal et al. 2006, 2007; Mori & Rich 2008; Kirihara, Miki & Mori 2014, 2017b), or in a major, gas-rich merger (e.g. D’Souza & Bell 2018; Hammer et al. 2018). This accreted event may also produce shell structures in the east and west of the M31 disc, and these shell structures have indeed been shown to be associated with GSS by spectroscopic observations (Escala et al. 2021; Dey et al. 2023). Comparing the observations of the M31 halo with those of the Galactic halo, it is suggested that disc galaxies such as the Galaxy and M31 have undergone different accretion and evolutionary paths, suggesting that individual stellar haloes differ in their star formation history and dwarf galaxy accretion times (e.g. Bullock & Johnston 2005; Cooper et al. 2010; Merritt et al. 2016; Monachesi et al. 2016; Harmsen et al. 2017, 2023). Thus, detailed observations of the M31 stellar halo are of great importance in gaining new insights into formation scenarios for the formation of galactic haloes.

However, since M31 is located at low galactic latitude, many foreground stars of the Galaxy make it difficult to distinguish M31 halo stars from such contaminations. Therefore, it is difficult to observe low surface brightness structures of the M31 halo. These Galactic foreground stars are expected to increase exponentially with decreasing galactic latitude, accounting for up to 70 per cent of all stars per deg|$^2$| (Martin et al. 2013). Therefore, although substructures have been discovered in M31, little progress has been made in estimating the distributions of their distances and line-of-sight velocities from observations and in simulating the origin of these substructures, except for the GSS and the stellar stream in the north-western part of the M31 halo (NW Stream; McConnachie et al. 2009).

To overcome this situation, a narrow-band filter (NB515) was developed for the Subaru Telescope/Hyper Suprime-Cam (HSC). With NB515, red giant branch (RGB) stars in the M31 halo can be efficiently separated from foreground dwarf stars in the MW based on differences in surface gravity between these stellar populations. In fact, for the NW Stream located at low galactic latitude (⁠|$\mathrm{b} \sim -15^{\circ }$|⁠), observations with NB515 succeeded to extract only M31 halo stars and estimate the distance to the stream (Komiyama et al. 2018b). From this result, it has become possible to constrain the accretion trajectory that reproduces the NW Stream (Kirihara, Miki & Mori 2017b). Using DDO51, an intermediate filter for the Mosaic Camera on the Kitt peak NationalObservatory 4 m Mayall telescope, RGB stars in the periphery of M31 were extracted, followed by spectroscopic observations using Keck/DEIMOS (e.g. Gilbert et al. 2006). These observations detected RGB stars at distances up to 175 kpc in the M31 halo (Gilbert et al. 2012). Therefore, the use of narrow-band filter such as NB515 for Subaru/HSC is powerful for understanding the physical properties of M31 substructures.

In this paper, we report observations and analyses of the M31 halo using the Subaru/HSC NB515. The HSC is a prime focus camera on the Subaru Telescope with a field of view diameter of 1.5 deg (Miyazaki et al. 2012, 2018; Furusawa et al. 2018; Komiyama et al. 2018a). We focus on the eight substructures already discovered in M31 and the properties of the halo from the data obtained with NB515. From this analysis, we demonstrate the effectiveness of NB515 by comparing the results with and without NB515, and also show the powerfulness of NB515 for future HSC analysis into deep broad-band data. The paper is organized as follows. In Section 2, we describe our observational and archival data and the method for data analysis of the M31 halo. The results of the distance estimation of the M31 substructure, the metallicity distribution, and the surface brightness profile of the halo are given in Section 3, and the discussion for these results is presented in Section 4. Finally, Section 5 concludes this paper.

2 DATA AND METHOD

2.1 Data

We use Pan-Andromeda Archaeological Survey (PAndAS; McConnachie et al. 2009, 2018) g- and i-band data combined with Subaru/HSC NB515 data as explained below to investigate the nature of the M31 stellar halo. The PAndAS data provide the wide and deep observational data for the M31/M33 stellar halo, so many studies have been conducted using the PAndAS data to understand the structure of the M31 stellar halo and to identify dwarf galaxies and stellar substructures in the M31 halo (e.g. McConnachie et al. 2009; Martin et al. 2013; Ibata et al. 2014). However, most of these studies assumed the statistical model to eliminate foreground contaminations (Martin et al. 2013). In this study, to reduce the uncertainties associated with such a model we use the NB515 filter for Subaru/HSC in combination with these PAndAS/g- and i-band data to select the M31 halo stars (see Sections 2.4 and 2.5) and to achieve the deep and homogeneous analysis for the M31 halo without contaminations.

The M31 halo was observed in 2015 and 2019 using the Subaru/HSC. Fig. 1 shows the observed regions in this study, where red and blue circles are the fields observed in 2019 (28 fields) and 2015 (5 fields), respectively (Komiyama et al. 2018b). These observed fields occupy about 50 deg|$^2$| of the M31 halo, and each field was observed multiple times to avoid CCD gaps and cosmic ray effects. The integration time varies from 960 to 1920 s, with most fields having integration times of 960 s, and its seeing is |$0{_{.}^{\prime\prime}} 50-1{_{.}^{\prime\prime}} 1$|⁠, with a mean value of |$0{_{.}^{\prime\prime}} 8$| and a standard deviation of |$0{_{.}^{\prime\prime}} 2$|⁠. Details of HSC observations are shown in Table 1.

The area that the PAndAS data covers, which can be accessed from the Canadian Astronomy Data Centre (CADC)¹ is shown in Fig. 1 within the outer black frame. PAndAS covers over 400 deg|$^2$| of the M31 halo and M33 as well, using the g and i bands, respectively. The integration time is 500 and 1160 s for g and i bands, and the data are available from the tip of the red giant branch (TRGB) of the M31 population to at least 3 mag fainter than TRGB. In the HSC observations, we performed dithering to compensate for the large and small gaps between the CCDs. Similar to the HSC observations, PAndAS also performed dithering by shifting the field of view of each shot by approximately |$10\,\mathrm{ arcsec}$|⁠. Although this could compensate for the small gap between each CCD, it was not enough to compensate for the large gap between the CCDs. Therefore, geometric gaps can be seen in the final spatial distribution (e.g. Fig. 7) in this study.

2.2 Reduction and photometry

The data observed by HSC are processed using the hscPipe 6.7 (Bosch et al. 2018). The hscPipe is based on the Large Synoptic Survey Telescope (Ivezic et al. 2008) pipeline and optimized for HSC. In the hscPipe, data reduction is performed on individual CCD, including de-bias, dark subtraction, flat-fielding, sky subtraction, and cosmic ray removal. After the reduction of individual CCD, the coordinates and flux scales of each CCD are calibrated, each frame is co-added, and finally, the source detection and photometry are performed to construct the source catalogue. In the hscPipe, it uses Pan-STARRS1 (PS1; Schlafly et al. 2012; Tonry et al. 2012; Magnier et al. 2013; Flewelling et al. 2020) for the photometric and astrometric calibration.

The final multicolour catalogue used in this study is constructed by cross-matching the HSC/NB515 catalogue with the PAndAS catalogue under the condition that the position of each star is |$\lt 0{_{.}^{\prime\prime}} 5$|⁠. In this study, we use only the sources that are determined to be stellar-like in both the hscPipe and PAndAS catalogues. In other words, we consider as stellar sources that are determined to have |$\mathrm{\tt {extendedness}} = 0$|⁠, which is indicators of successful star–galaxy classification at depths from |${\it i}\sim 23~{\rm mag}$| (Aihara et al. 2018) in hscPipe and have a probability of being a point source within |$1\sigma$| or |$2\sigma$| in both g- and i-bands of the PAndAS catalogue.

Based on the dust maps of Schlegel, Finkbeiner & Davis (1998), we apply Galactic extinction correction for our catalogue. The extinction-corrected magnitudes for each band are

where the subscript zero means the extinction-corrected magnitude. The coefficient of the extinction-corrected magnitude of NB515 is obtained by multiplying the interstellar absorption curve of Fitzpatrick (1999) with |$R=+3.1$| by the spectral energy distribution (SED) of a G-type star (⁠|$\mathrm{\mathit{ T}_{eff}}=7000$| K, |$\log {\mathrm{Z}}=-1$|⁠, |$\log {g}=4.5$|⁠) and integrating it using the NB515 response curve. For |$g_0$| and |$i_0$|⁠, we use the coefficient from PAndAS studies (e.g. Ibata et al. 2014; McConnachie et al. 2018).

The detection completeness varies from region to region depending on the signal-to-noise ratio (S/N) and the crowding of objects. We evaluated the detection completeness of NB515 by embedding and detecting artificial stars in the images. To perform this test, we developed a python module injectStar.py.² Details of this module and the artificial star test are described in Appendix  A. Fig. 2 shows the apparent magnitudes corresponding to the 50 per cent completeness of NB515 for each |$0.1 \times 0.1~{\rm degree}^2$| region. The average magnitude of the 50 per cent completeness in the NB515 band is 23.21 mag. In this map, the region where seeing is slightly worse (the centre coordinate of the tangential plane are (3,1) for |${\rm PFS}\_{\rm FIELD}\_26$|⁠, |$(-1,-2)$| for |${\rm PFS}\_{\rm FIELD}\_31$|⁠, |$(0,-3)$| for |${\rm PFS}\_{\rm FIELD}\_32$|⁠), the 50 per cent detection completeness magnitude is the almost same as the other regions. Detection completeness for the PAndAS survey was calculated and modelled by Martin et al. (2016). They showed the completeness |$\eta (m)$| in equation (2).

where m is the apparent magnitude of the g and i bands, |$m_{50}$| is the apparent magnitude of 50 per cent completeness, A and |$\rho$| are the constant values. In the g band, |$m_{50} = 24.88$|⁠, |$A = 0.94$|⁠, and |$\rho = 0.65$|⁠, while in the i band, |$m_{50} = 23.88$|⁠, |$A = 0.93$|⁠, and |$\rho = 0.74$|⁠. In this study, the completeness model is assumed to be uniform in the overall region. By fitting the results of the artificial star test (the ratio of the number of detected artificial stars to the number of injected artificial stars) to equation (2), we obtain |$m_{50} = 23.21$|⁠, |$A = 0.99$|⁠, and |$\rho = 0.30$| for NB515.

2.3 Colour–magnitude diagram

Fig. 3 shows the de-reddened colour–magnitude diagram (CMD) shown as the Hess diagram in the entire region analysed in this work, and the Dartmouth isochrones (Dotter et al. 2008) with |$[\alpha /\mathrm{Fe}] = 0.0$|⁠, age of 13 Gyr, and different metallicities of |$[\mathrm{Fe/H}]=0.00, -0.70, -1.50, -2.00$| which are scaled to the distance of M31 (776; Dalcanton et al. 2012). The yellow dashed line shows the 50 per cent completeness, and the red crosses are the average photometric errors at |$(g-i)_0 = 1.0$|⁠. From this figure, we can identify various stellar populations. The RGB stars in the M31 halo, the target of this study, are distributed from |$((g-i)_0, i_0) \sim (2.0, 21)$| to |$\sim (0.5, 25)$|⁠. The region of |$18 \lessapprox i_0 \lessapprox 22$| at |$(g-i)_0 \sim 2.2$| is occupied by the main sequence (dwarf) stars in the foreground Galactic disc, while the region of |$18 \lessapprox i_0 \lessapprox 22$| at |$(g-i)_0 \sim 0.5$| shows the main sequence turn-off (MSTO) stars in the Galactic halo. These stellar populations indicate that the RGB stars in M31 overlap with the foreground Galactic stars in the CMD. Previous PAndAS works have used data with |$i_0\lt 23$| to account for the detection completeness and photometric errors (e.g. Martin et al. 2013; Bate et al. 2014), so we also use the stars with |$i_0\lt 23$| for the analysis.

2.4 NB515 and colour–colour diagram

Fig. 3 also shows that the distribution of the M31 RGB stars is overlapping with that of the Galactic dwarf stars in the CMD, making it difficult to distinguish these stellar populations from only the |$(g_0, i_0)$| information. For this purpose, we use the Subaru/HSC narrow-band filter NB515, which is sensitive to the surface gravity of a star: the intensity of spectral absorption in the MgH band at 521 nm and in the Mgb triplet at 518 nm strongly depend on the stellar surface gravity (Öhman 1934; Thackeray 1939). Based on this property, Majewski et al. (2000) proposed a method to separate RGB stars from dwarf stars using DDO51 which is an intermediate-band filter with a central wavelength of 515 nm and bandwidth (FWHM) of 12.3 nm (Clark & McClure 1979). Since the MgH + Mgb absorptions are covered by DDO51, there is a difference in the DDO51 magnitude between the RGB and dwarf owing to different surface gravities: the colour index |$({\it M}-{\it DDO51})$|⁠, where M is a broad-band magnitude in the Washington system, provides an index of the surface gravity (Majewski et al. 2000).

Following this property of DDO51, NB515 for HSC is made with a central wavelength of 515 nm and an FWHM of 7.7 nm by our group (PI: M. Chiba). In Fig. 4, the solid black line shows the sensitivity curve of NB515 only, and the dotted black line shows the sensitivity curve of NB515 multiplied by the atmospheric transmittance, the reflectance of the primary mirror, the transmittance of the prime focus optics, the quantum efficiency of the CCD, and the transmittance of the dewar window. Fig. 4 also shows stellar model spectra with |$\mathrm{\mathit{ T}_{eff}} = 5000$| K, |$[\mathrm{C/M}] = 0.00$|⁠, and |$[\alpha /\mathrm{M}] = 0.25$| BOSZ stellar model spectra³ (Bohlin et al. 2017). The red and green lines show dwarf stars with |$[\mathrm{Fe/H}]$| of 0.00 and |$\log {g}$| of 4 and 5, respectively, while the orange and blue lines show giant stars with |$[\mathrm{Fe/H}]$| of 0.00 and |$-2.00$|⁠, respectively, and |$\log {g} = 1$|⁠. In Fig. 4, the absorption lines that overlap with the sensitivity curve of NB515 correspond to MgH + Mgb, and it can be confirmed that these lines for RGB stars (orange and blue) are weaker than those for dwarf stars (red and green). Following the method proposed by Majewski et al. (2000), we use three photometric bands (g, i, NB515), which cover a similar wavelength range as those photometric bands used in Majewski et al. (2000), to extract M31 RGB stars on the colour–colour diagram.

Fig. 5 shows the colour–colour diagram of the entire analysis region and the theoretical curves calculated by BOSZ model spectra with the properties of the dwarf and RGB stars. The blue and orange solid lines are the RGB-like model spectra with |$[\mathrm{M/H}] = -2.00$|⁠, |$\log {g} = 1.0$|⁠, and with |$[\mathrm{M/H}] = 0.00$|⁠, |$\log {g} = 1.0$|⁠, respectively. The green and red dashed lines correspond to the dwarf-like model spectra with |$[\mathrm{M/H}] = 0.00$|⁠, |$\log {g} = 4.0$|⁠, and with |$[\mathrm{M/H}] = 0.00$|⁠, |$\log {g} = 5.0$|⁠, respectively. It is clearly seen that the distribution of RGB stars (enclosed by a solid cyan rectangle) is different from that of dwarf stars (enclosed by a dashed cyan polygon). Note that the metal-rich giant (orange solid line in Fig. 5) intersects the distribution of many dwarf stars (a cyan dashed polygon) with |$({\it g}-{\it i})_0 \gt 2$|⁠. Therefore, our NB515-based selection cannot work well on the metal-rich side, so to account for this, we cut these red stars in our analysis (metallicity profiles and surface brightness profiles in the halo of M31; Section 3).

2.5 Selection of the M31 RGB stars

In this subsection, we derive the dwarf probability (RGB probability) based on narrow-band information, |$p_{\mathrm{dwarf, {\it NB}}}$| (⁠|$p_{\mathrm{RGB, {\it NB}}}$|⁠) and dwarf probability (RGB probability) based on galactic latitude information |$p_{\mathrm{dwarf, {\it lat}}}$| (⁠|$p_{\mathrm{RGB, {\it lat}}}$|⁠) for each star. From the combination of these probabilities, we calculate the probability of being an RGB star in M31, |$p_{\mathrm{M31}}$|⁠, for each star.

First of all, we assume that the distance to disc of M31 is 776 kpc (Dalcanton et al. 2012) and all M31 RGB stars are located within a virial radius of |$\sim$|260 kpc (Seigar, Barth & Bullock 2008). Then, while the M31 RGB stars are distributed mostly at the fainter side of |$i_0 \gtrsim 21$| mag of the CMD, many foreground Galactic dwarf stars are available at the brighter side, in |$i_0 \lt 20$| mag, where there exist no M31 RGBs. The left panel of Fig. 6 shows the colour–colour diagram of stars with |$i_0 \lt 20$|⁠, which is regarded as the colour–colour diagram of only foreground dwarf stars. We bin the data between |$0.0 \le ({\it g}-{\it i})_0 \le 3.0$| in |$({\it g}-{\it i})_0$| bins of width 0.2 mag. In each bin, we compute the mean and standard deviation of |$({\it NB515}-{\it g})_0$|⁠. By interpolation of these mean values using the python/scipy package, we construct the function that represents the distribution of dwarf stars (Dwarf Ridge Line; |$f_{\mathrm{dwarf}}((g-i)_0)$|⁠). The left panel of Fig. 6 shows Dwarf Ridge Line |$f_{\mathrm{dwarf}}((g-i)_0)$| as a solid red line, where the red dots and error bars depict the mean values and 3 times standard deviations, respectively, in each bin.

To easily derive the |$p_{\mathrm{dwarf},{\it NB}}$|⁠, we introduce |$\Delta _{{\it NB}} \equiv ({\it NB515}-{\it g})_0-f_\mathrm{dwarf}(({\it g}-{\it i})_0)$|⁠, which makes the mean of the distribution of foreground dwarf stars a straight horizontal line on the diagram. The middle panel of Fig. 6 shows the |$({\it g}-{\it i})_0$|–|$\Delta _{{\it NB}}$| diagram, where the vertical axis of the left panel of Fig. 6 is replaced by |$\Delta _{{\it NB}}$|⁠. We assume that the probability distribution of the dwarf sequence at a given |$(g-i)_0$| is a Gaussian distribution with mean |$\mu = 0$|⁠, and standard deviation |$\sigma = \sigma _{({\it NB515}-{\it g})_0}$| in |$\Delta _{{\it NB}}$| direction, where the standard deviation |$\sigma _{({\it NB515}-{\it g})_0}$| is calculated for each bin of 0.2 mag. Note that the photometric errors of dwarf stars (⁠|${\it i}_0\lt 20$|⁠) are sufficiently small (mean value is |$\lt 0.01$|⁠) that |$\sigma _{({\it NB515}-{\it g})_0}$| does not include photometric error.

Performing the linear interpolation of these normal distributions, we construct the probability distribution |$P_{\mathrm{dwarf}}$| in which individual stars are likely to be the dwarf star on the |$({\it g}-{\it i})_0$|–|$\Delta _{{\it NB}}$| diagram. The right panel of Fig. 6 shows this distribution.

Based on this right panel and the photometric errors of individual stars, we calculate the dwarf probability |$p_{\mathrm{dwarf},{\it NB}}$| for each star. First, we assume that each star is distributed as a two-dimensional Gaussian distribution based on the photometric errors in colour |$\sigma _{(g-i)_0}$| and |$\sigma _{({\it NB515}-{\it g})}$|⁠. Then, the distribution function |$\mathcal {N}(({\it g}-{\it i})_0,\Delta _{{\it NB}})|(g-i)_{0,n}, \Delta _{{\it NB },n})$| of the n-th star with photometric errors of |$(\sigma _{({\it g}-{\it i})_{0,n}}, \sigma _{\Delta _{{\it NB},n}})$| at |$(({\it g}-{\it i})_{0,n}, \Delta _{{\it NB},n})$| in the colour–colour diagram is expressed as

where

and |${\rm Cov}((g-i)_0,\Delta _{NB})$| is the covariance between |$(g-i)_0$| and |$\Delta _{NB}$|⁠. It is noted that when |${\rm Cov}((g-i)_0,\Delta _{NB})$| is larger than |$\sigma _{(g-i)_{0,n}}^2$| or |$\sigma _{\Delta _{{\it NB},n}}^2$|⁠, V is no longer a positive semidefinite matrix, and we cannot calculate the bivariate normal distribution (equation 3). Therefore, when |${\rm Cov}((g-i)_0,\Delta _{NB})$| is greater than |$\sigma _{(g-i)_{0,n}}^2$| or |$\sigma _{\Delta _{{\it NB},n}}^2$|⁠, |${\rm Cov}((g-i)_0,\Delta _{NB})$| is assumed to be 0, because when the photometric error is small, it is expected that the distribution of individual stars on the colour–colour diagram will not differ significantly.

The dwarf probability |$p_{\mathrm{dwarf},{\it NB}}$| for each star is calculated by integrating the product of the probability distribution of dwarf stars |$P_{\mathrm{dwarf}}$| shown in the right panel of Fig. 6 and the distribution of individual stars shown in equation (3) as follows:

and normalizing these calculated values. The RGB probability |$p_{\mathrm{RGB},{\it NB}}$| of each star is calculated to be |$p_{\mathrm{RGB},{\it NB}} = 1-p_{\mathrm{dwarf},{\it NB}}$|⁠, assuming that the stars distributed in the colour–colour diagram are only Galactic dwarf stars and M31 RGB stars.

From the colour–colour diagram (see Figs 5 and 6), it is evident that dwarf stars generally have |$\Delta _{{\it NB}}\gt 0$|⁠, while RGB stars have |$\Delta _{{\it NB}} \lt 0$|⁠. Therefore, we take into account this effect by multiplying |$p_{\mathrm{RGB,{\it NB}}}$| with the step function |$f_{\mathrm{step}}$|⁠,

to make the value of |$p_{\mathrm{RGB},{\it NB}}$| smaller when |$\Delta _{{\it NB}} \gt 0$|⁠.

The number of foreground stars is expected to change over our |$\sim 50$| deg|$^2$| field of view as a function of latitude (see e.g. Tanaka et al. 2010; Martin et al. 2013). We incorporate this spatial effect into our probabilities by computing |$p_{\mathrm{RGB},lat}$| as follows. First, in the same way as Tanaka et al. (2010), we fit the spatial distribution of foreground-like stars with |$i_0\lt 20$| with a function of galactic latitude, b,

where A, B, and C are the constant values. In this function, the returned value (dwarf probability of galactic latitude) is 1 at |${\it b}=0$|⁠. Then, based on the spatial information for each star, we calculate the dwarf probability |$p_{\mathrm{dwarf},{\it lat}}$| of each star. The RGB probability |$p_{\mathrm{RGB},{\it lat}}$| based on the latitude is calculated as |$p_{\mathrm{RGB}, lat} = 1-p_{\mathrm{dwarf}, lat}$| under the assumption that there are only RGB stars and dwarf stars as stellar populations in this analysis. Finally, based on the calculated RGB probabilities of each star based on its NB515 and galactic latitude, the probability of being an RGB star |$p_{\mathrm{M31}}$| is calculated according to the following equation:

where the denominator is the normalization constant.

Fig. 7 shows the spatial distribution (top panels) and CMD (bottom panels) for all stars (leftmost panels) and those with M31 RGB probability |$p_{\mathrm{M31}} \gt 0.9$| (inner-left panels). It is clear that the CMD for the stars with high |$p_{\mathrm{M31}}$| show a more prominent RGB sequence than those without the selection by |$p_{\mathrm{M31}}$|⁠. In the spatial distribution, although the number of objects decreases when the selection by |$p_{\mathrm{M31}}$| is applied, we can confirm the already discovered substructures in the M31 halo (e.g. GSS, Stream C, and Stream D; Ibata et al. 2001, 2007; Ferguson et al. 2002). Therefore, by using the RGB probability |$p_{\mathrm{M31}}$|⁠, the clear structure of the M31 halo is made available without the influence of contamination from the Galactic foreground dwarf stars.

Using |$p_{\mathrm{M31}}$|⁠, it is possible to extract a large number of M31 RGB stars. The inner-right bottom panel in Fig. 7 shows the CMD for the stars with high |$p_{\mathrm{M31}}$| and low galactic latitude (⁠|$\eta \gt 2.5$| on M31-centric coordinates). In this panel, we confirm that there are stars at |$(({\it g}-{\it i})_0) \sim (2.5, 22)$|⁠, which look like remaining Galactic dwarf-like stars. As seen in Fig. 5, it is difficult to separate the dwarf stars from the RGB stars at |$(g-i)\gt 2.5$|⁠. This causes some (red) dwarf stars to have high |$p_{\mathrm{RGB}}$|⁠. It should be noted, however, that the contrast of the inner-right bottom panel in Fig. 7 is different from the other CMDs, so we expect that the number of the remaining dwarf stars is small. In Section 4.1, we verify the accuracy of the M31 RGB probability |$p_{\mathrm{M31}}$| calculated in this study, confirm that the accuracy of |$p_{\mathrm{M31}}$| is |$\sim 90~{{\ \rm per\ cent}}$|⁠.

In order to remove these remaining dwarf stars with high |$p_{\mathrm{M31}}$| and investigate the pure properties of the M31 stellar halo, we further set an RGB box to select the most likely RGBs on the CMD based on isochrones, where the age is 13 Gyr, |$[\alpha /\mathrm{Fe}] = 0.0, Y = 0.245 + 1.5Z$|⁠, and |$[\mathrm{Fe/H}] = -2.5$| to |$+0.5$|⁠. To achieve this, we assume isochrones to be distributed at a distance of |$785\pm 100$| kpc, assuming the spread of the M31 stellar halo. The fainter boundary of |$i_0$| of the RGB box is set to |$i_0 = 23$| by considering the completeness. The reddest of |$(g-i)_0$| is set to |$(g-i)_0 = 2.5$|⁠, considering that completeness decreases on the red side and that RGB stars overlap dwarf stars on the red side. The CMD in the rightmost bottom panel of Fig. 7 shows the RGB box as the red box. Hereafter, the objects with |$p_{\mathrm{RGB}} \gt 0.9$| in this RGB box are called NRGB stars (basically NB515-based selected RGB stars). The NRGB stars are used to construct the photometric metallicity distribution described in Section 3.2.1 and the surface brightness profile described in Section 3.2.2.

2.6 Algorithm for distance estimation

In metal-poor stellar systems in which individual stars can be resolve, the brightness of the TRGB can be used as a distance indicator. Starting with Lee (1993), the absolute i-band magnitude of TRGB was found to vary by only 0.1 mag in stellar systems with |$[\mathrm{Fe/H}] \lt -0.7$| (e.g. Lee 1993; Bellazzini, Gennari & Ferraro 2005; Rizzi et al. 2007). This feature leads to a sharp drop at magnitudes brighter than the apparent magnitude of TRGB for stellar systems with low metallicity and a single stellar population. The luminosity function of such RGB stars with apparent magnitudes, m, can be modelled as

where a and b are constants. The value of parameter a should originally be estimated using the initial mass function, but since it has been treated as a free parameter in RGB distance estimation (e.g. Makarov et al. 2006; Conn et al. 2011; Tollerud et al. 2016), so it is regarded as a free parameter in this study as well.

Such distance estimation methods are effective for stellar systems with low metallicity and a single population. However, for stellar systems with metal-rich and/or a wide range of metallicity, the TRGB brightness is not constant and no clear step-like distribution appears in the luminosity function, even for systems with a sufficient number of objects. Therefore, for stellar systems with metal-rich and/or a wide metallicity dispersion, it is not possible to fit with equation (8). To overcome this problem, a distance estimation method is developed by constructing a model CMD based on the RGB isochrones on the CMD and fitting the model to the observed data (Conn et al. 2016). Some substructures that are analysed in this study include stellar systems that are known to have metal-rich and a wide metallicity dispersion. Therefore, the method developed by Conn et al. (2016) was optimized and applied to the data in this analysis to estimate the distance. In this section, we briefly introduce the method used in Conn et al. (2016) to construct the model CMD and the distance estimation method, as well as the optimization of the method to our analysed data.

Following the procedure of Conn et al. (2016), we perform the model fitting and distance estimation, applying likelihood calculations and the Markov chain Monte Carlo method (MCMC) to the observed data. To construct the colour–magnitude model, we prepare the isochrones with the metallicity of |$-2.5 \le [\mathrm{Fe/H}] \le 0.5$| and the age of 10 Gyr. The helium mass ratio is assumed to be |$Y = 0.245 + 1.5Z$|⁠, and the |$\alpha$|-element ratio is assumed to be |$[\alpha /\mathrm{Fe}] = 0.0$|⁠. In constructing the model, the distribution of the number of RGB stars on each isochrone is obtained and then the distribution is interpolated. The distribution of the number of RGB stars on each isochrone is constructed using the above equation (8) and the following equation (9). First, the number of RGB stars at a given magnitude m is assumed to vary according to equation (8). In this equation, |$m_{\mathrm{TRGB}}$|⁠, a, and b are the parameters of the model. |$m_{\mathrm{TRGB}}$| is the magnitude of TRGB, a is the slope of the luminosity function, b is the model of the foreground/background objects. The parameter b is a constant, and m is the apparent magnitude in the isochrone. In the Dartmouth isochrones, the available isochrone data set is given by the absolute magnitude |$M_{\mathrm{RGB}}$|⁠. Therefore, the absolute magnitude |$M_{\mathrm{RGB}}$| needs to be converted to apparent magnitude |$m_{\mathrm{RGB}}$|⁠. Here, we use the distance to the stellar system in kpc, d, and d is a parameter to be estimated during model fitting. In the observed stellar system, the metallicity distribution is not uniform metallicity distribution. Therefore, it is necessary to reflect the gradient of the number of objects according to the metallicity distribution in the model CMD of RGB stars. The contribution of an isochrone with a certain metallicity |$[\mathrm{Fe/H}]_{\mathrm{iso}}$| to the model CMD is given by

assuming that the metallicity distribution of the stellar system is the Gaussian distribution. The parameters |$[\mathrm{Fe/H}]_0$| and |$\omega _{\mathrm{RGB}}$| are the mean metallicity and standard deviation of the metallicity of the stellar system. Multiplying equations (8) and (9) determines the distribution of the number of RGB stars on each isochrone. Therefore, the change in the number of RGB stars in an isochrone with a certain metallicity |$[\mathrm{Fe/H}]_{\mathrm{iso}}$| becomes

A model CMD is constructed by three-dimensional linear interpolation of equation (10). Fig. 8 shows the model around the centre of M31 in the GSS created by applying the distance estimation method. The red dots are the observed data, and it can be confirmed that the model (black shaded) reproduces the observed RGB stars (the distribution from |$((g-i)_0 \sim 2, i_0 \sim 21)$| to |$((g-i)_0 \sim 1.5, i_0 \sim 22))$|⁠. Conn et al. (2016) use the model constructed above to estimate distances by calculating the likelihood |$\mathcal {L}_{\mathrm{CMD},n}$| of n-th observed data and applying it to MCMC.

In addition to the above, this study considers the detection completeness and RGB probability |$p_{\mathrm{M31}}$| for the calculation of the likelihood. For the n-th star, we calculate the detection completeness (⁠|$\eta ({\it g}_{0,n}), \eta ({\it i}_{0,n}), \eta ({\it NB515}_{0,n})$|⁠) from equation (2) and the probability |$p_{\mathrm{M31},n}$| (see Sectioon 2.5). We assume that the contribution of the n-th data to the CMD of the substructure is weighted by |$p_{\mathrm{M31}}$| and the inverse of the detection completeness. Under this assumption, we can derive the likelihood by directly multiplying these values to the likelihood |$\mathcal {L}_{\mathrm{CMD},n}$| from the model CMD, so the likelihood is calculated as follows:

Assuming that the prior distribution of the parameters is the uniform distribution of the intervals shown in Table 2, the python module emcee (Foreman-Mackey et al. 2013) is used to estimate the distance for each subregion of each substructure. We initialize the sampler with 10 walkers, with 100 000 steps after a burn-in period of 400 000 steps to ensure a good sampling of the posterior distributions.

2.7 Metallicity estimation

The photometric metallicities of individual stars are estimated by comparing the positions of stars on the CMD with the Dartmouth isochrones (Dotter et al. 2008) with given metallicities as follows. To estimate the photometric metallicity for each star, we construct 31 isochrones with |$-2.5 \lt [\mathrm{Fe/H}] \lt +0.5$| in 0.1 dex intervals. To simplify the comparison with the previous photometric metallicity studies (e.g. Ibata et al. 2014; McConnachie et al. 2018), we assume that the age is 13 Gyr and alpha-abundance is |$[\mathrm{\alpha /Fe}]=+0.0$| which is the same as the previous studies. It is noted that the alpha abundance of the M31 halo is not zero (⁠|${\rm [\alpha /Fe]} \gt 0.0$| Wojno et al. 2023), so our metallicity estimation may underestimate by |$\lt 0.2$| dex, systematically.

For the metallicity estimation, we also assume the distance to M31 to be 776 kpc (Dalcanton et al. 2012) and construct the metallicity model of RGB stars on the CMD by using the radial basis function (Rbf) interpolation of the python/scipy package. Fig. 9 shows the interpolated RGB metallicity model. The metallicity of each star at a given position |$((g-i)_0, i_0)$| on the CMD is estimated by interpolating [Fe/H] at that position in the metallicity model. By this method, we estimate the photometric metallicity of all point sources.

3 RESULTS

In this section, we describe the results of our analysis of the M31 halo, based on the catalogue constructed in the previous section. We describe the properties of the substructures in the M31 halo in Section 3.1, and the natures of the stellar halo in the entire analysis region in Section 3.2.

3.1 The properties of M31’s substructures

3.1.1 The spatial distribution for M31’s substructures

The spatial distribution of the NRGB stars is shown in Figs 10 and 11. Fig. 10 shows the spatial distribution of stars with high |$p_{\mathrm{M31}}$| (i.e. including those red stars excluded in the NRGB colour–magnitude selection) and the substructures that have already been found. Fields of substructures are indicated by enclosing them with polygons which are determined visually. Fig. 11 shows the spatial distribution smoothed with a Gaussian kernel over 3 pixels (⁠|$= 1.8$| arcmin) and divided by [Fe/H] into four ranges (−0.6 to −0.1, −1.1 to −0.6, −1.6 to −1.1, and −2.1 to −1.6). In these maps, we can see the substructures (e.g. GSS, Eastern and Western Shell, Stream C and D, and NW Stream) that have already been found in previous studies. In Section 3.1.2, we estimate the parameters (e.g. distance and metallicity) of the substructures that have already been found.

The substructures already found show a spatial distribution in good agreement with previous studies. In addition, we can identify the two overdense structures on the north-west halo at |$(\xi , \eta) \sim (-1^{\circ }, 2{_{.}^{\circ}} 5)$| in the most metal-rich map and at |$(\xi , \eta) \sim (-3^{\circ }, 1^{\circ })$| in the most metal-poor map (hereafter, we call them the Metal-Rich Cloud and the Metal-Poor Cloud). Moreover, we can confirm one stream-like structure on the south-east halo at |$(\xi , \eta) \sim (1{_{.}^{\circ}} 5, -1^{\circ })$| (hereafter, we call this stream-like structure the SE Stream) from the visual inspection of the Gaussian kernel-smoothed maps (the red boxes in Fig. 11 indicates the newly discovered substructures). The blue boxes in Fig. 11 are assumed to be off-substructure fields to calculate the S/N of newly discovered substructures. Fig. 12 is a zoomed-in view of the newly discovered substructure field. The red and blue boxes are the same as Fig. 11. The contours represent an S/N of 9 (dark red), 6 (light green), and 3 (navy), where S/N is estimated through a formula |$=$|(smoothed map—mean of off-substructure field)/standard deviation of off substructure field (Yang et al. 2023). As can be seen from this figure, the new substructures (red boxes) have a S/N of 6 or higher for the off-substructure fields (blue boxes), and these high S/N are the almost same as Stream C (⁠|$(\xi ,\eta) \sim (2{_{.}^{\circ}} 5,-1^{\circ })$| in the middle panel of Fig. 12) and NW Stream 1 (⁠|$(\xi ,\eta) \sim (-0{_{.}^{\circ}} 5,4^{\circ })$| in the right panel of Fig. 12).

Fig. 13 shows the Hess diagram with normalized number of stellar objects of Metal-Rich, Metal-Poor Clouds (fields defined by the red boxes in Fig. 11), off Metal-Rich/Metal-Poor Clouds region (defined by the blue boxes in Fig. 11), Western Shell, and NW Stream. This figure shows the well-populated RGB sequence for Metal-Rich/Metal-Poor Clouds, in contrast to the off-clouds. The spatial distribution of both clouds slightly overlaps each other and it seems that most part of these new overdensities are an outward extension of the Western Shell. Also, for the Metal-Rich Cloud, NGC205-loop and North Spur are nearby, but this new one is located outside of these known structures. Only Andromeda XVI overlaps the Metal-Rich Cloud, but this dwarf galaxy is more metal-poor (⁠|$\mathrm{[Fe/H]} \sim -2$|⁠; Letarte et al. 2009; Collins et al. 2014, 2015) than the Metal-Rich Cloud (⁠|$-0.6\lt [\mathrm{Fe/H}]\lt -0.1$|⁠). Therefore, although it is difficult to distinguish this cloud from Andromeda XVI using only spatial distribution, metallicity (i.e. the colour–magnitude information) can separate them. In addition to the Metal-Rich Cloud, the photometric metallicity distribution of the Western Shell seems to be wider than that of Metal-Poor Cloud in Fig. 13, so we can consider that this substructure is not part of the Western Shell, according to the metallicity information. Based on the above, we consider these substructures to be new discovered ones. In Section 4.3, we describe distance estimation results for these clouds, and discuss the origin of these structures.

Fig. 14 shows the Hess diagram with the normalized number of stellar objects of the SE Stream, off-SE Stream, Stream C, and Stream D. In this figure, the RGB sequence of the SE Stream extends to the red side, and metal-poor sequences such as Stream C and D also can be seen. In particular, the CMD of the off-Stream has no noticeable sequence within the area enclosed by the red dashed lines, while the SE Stream has a well-populated RGB sequence even within the red dashed lines. Therefore, it can be confirmed that there are stellar populations that are a combination of off-Stream components (i.e. field halo stars) and Stream C and/or D components (i.e. accreted stars; shown in the left panels of Fig. 14). The South-Eastern Shell (SE shell; Gilbert et al. 2007) is very close to this stream. However, since this structure is slightly outside of the SE shell and is more metal-poor (⁠|$[\mathrm{Fe/H}]\sim -1.3$|⁠) than the SE shell (⁠|$[\mathrm{Fe/H}]\sim -0.50$|⁠), we can consider that the SE Stream is different from the SE shell. In Section 4.3, we describe distance estimation results for these clouds, and discuss the origin of this structure, as well as Metal-Rich/Metal-Poor Clouds.

3.1.2 Line-of-sight distances to M31’s substructures

We apply the method described in Section 2.6 to estimate the distances of the halo substructures of M31. Fig. 10 shows the density map of the NRGB stars. The solid colour-coded lines are the outer frame for each of the M31 substructure in the analysis region. The dashed colour-coded lines correspond to the subregions into which the substructures are divided for the distance estimation: to derive the distance and its gradient in each substructure, we perform distance estimation for each subregion. Note that an overlapping system of fields was implemented such that a given subregion XY (e.g. GSS12 shown in Table 3) contains the stars from both subregion X and subregion Y. In addition, it should be noted that all stars with |$0.5 \lt ({\it g}-{\it i})_0\lt 2.5$| and |$20\lt {\it i}_0\lt 22$| are used in this analysis when the distance estimation. Finally, we should note that the regions for substructures are determined visually, so they may deviate slightly from the regions estimated in previous studies.

Table 3 is the results of the estimation for each of the M31 substructures, showing the distance, metallicity, and metallicity deviation, and their random errors (corresponds to 68 per cent Bayesian credible interval for posterior distribution). In addition to random error, we estimate the systematic error of the distance as follows. As described in Section 2.6, the distance estimation is performed by model fitting on the CMD, in which there exist several systematic errors. First, according to Tonry et al. (2012), Pan-STARRS 1 has a systematic error of 0.02 mag (⁠|$\sim 7$| kpc) for each g and i bands. PAndAS data used the PanSTARRS photometry to correct for zero-point uncertainties in the survey region, so these systematic photometric errors are included as total systematic errors of distance estimation. Second, as described in Section 2.1, we correct for the magnitude based on the dust map by Schlegel et al. (1998), and extinction correction is considered to have a systematic error of about |$10~{{\ \rm per\ cent}}$|⁠, which corresponds to a magnitude of 0.015 mag (⁠|$\sim$|5 kpc). Third, in this distance estimation method using the Dartmouth isochrones, we find that the distance changes by a maximum of 5 kpc, if we change the age assumption to 2 Gyr (⁠|$10 \pm 2$| Gyr). Finally, Serenelli et al. (2017) showed that RGB stars in such a theoretical curve contain a systematic error of 0.05 mag (⁠|$\sim 10$| kpc). Adding these four factors, we estimate the total systematic error in the estimated distance to be |$\sim 27$| kpc.

Fig. 15 shows the results of distance estimation for GSS. The leftmost panel is the spatial distribution with points where these are centres of stars within the distance-estimated region weighted by |$p_{M31}$|⁠. The inner-left panel is the result of the distance as a function of |$\eta$| coordinate. The inner-right is the estimated metallicity of stellar systems, and the rightmost panel shows the standard deviation of the metallicity. In each panel, red dots show the results of this study, and green and blue dots correspond to the results of Conn et al. (2016) and McConnachie et al. (2003), respectively. In the three figures to the right, the results of this study as red points are displayed by violin plots, where the shaded region represents the shape of the posterior distribution within the 68 per cent Bayesian credible interval. A black dot shows systematic error in the estimation of this study. The results for GSS are generally consistent with previous studies, and we can confirm the same distance gradient that is reported by previous studies. In addition, the results of estimated metallicity and its standard deviation are also almost consistent with those of Conn et al. (2016). While the overall trend is generally consistent, some subregions (e.g. GSS2, GSS3, GSS4) are estimated to be more metal-rich than those of Conn et al. (2016). This discrepancy may be attributed to the limited number of objects used in the analysis. However, when using objects from both subregions, such as GSS23 and GSS34, which increases the sample size, the results show a good agreement with those of Conn et al. (2016). While generally consistent with the metal amounts estimated in previous works, when compared to recent spectroscopic results (e.g. Escala et al. 2021) it shows that the metallicity distribution of GSS is more metal-poor (⁠|$[{\rm Fe/H}] \lt -0.8$|⁠) than our result as one moves away from the centre of M31. This result may be due to the lack of correction for the kinematically hot stellar halo of M31 in the analysis of photometric data. Therefore, it is necessary to confirm these metallicity peaks using future follow-up spectroscopic observations as Subaru/Prime-Focus Spectrograph (PFS; Takada et al. 2014).

Fig. 16 shows the results of the distance estimation for the Stream C and the Stream D, similar to Fig. 15. The red points in each figure show the estimation results and errors of this study, and the green dots indicate the results of Conn et al. (2016). Although our results are basically in agreement with previous results, we identify the presence of significant distance gradients for both streams clearly. We think that these trends can be more clearly seen thanks to the addition of the narrow-band data in our selection of candidate member stars. In Conn et al. (2016), stars other than the M31 substructure are treated statistically based on the method of Martin et al. (2013), in which the contamination model was constructed based on regions where the M31 substructure was undetected. In contrast, we assign each star to be a contamination source or not, based on its NB515 magnitude (see Section 2.5). Therefore, the difference in treating contamination sources is likely to be reflected in the distance estimation results and their precisions.

Fig. 17 shows the distance estimation results for the NW Stream. The red points and green dots are the results of this study and Komiyama et al. (2018b), respectively. In the region at |$\eta \gt 3$|⁠, the results agree with those of Komiyama et al. (2018b) within the error range. However, Komiyama et al. (2018b) concluded that the NW Stream is fully located behind M31 from their distance estimation using Red Clump (RC) stars. On the other hand, our analysis suggests that the NW Stream is distributed in front of M31 at |$\eta \lt 3$| (where the distance is below 776 kpc). We think this difference occurs due to the choice of different distance indicators and whether contamination removal is performed. Komiyama et al. (2018b) performed the estimation using RC, which differs from this study with RGB stars. In addition, they did not use NB515, because the NB515 data are shallower than the magnitude of RC in the M31 halo. RC (⁠|$\delta i_0 \sim \pm 0.05$|⁠) used in Komiyama et al. (2018b) is a larger intrinsic (⁠|$\sim$| distance) error than that of RGB stars (⁠|$\delta i_0 \sim \pm 0.015$| mag), and they do not take into account the contamination sources (e.g. foreground stars and background unresolved galaxies). Considering that there is |$\sim$|20 kpc of systematic errors in addition to the random errors, the total errors in distance estimation are large and it is not possible to rule out the possibility that the stream is in front of M31. Also, the RC method by Komiyama et al. (2018b) contains the effect of contamination sources, despite the fact that foreground stars are severely contaminated in the region where the NW Stream is located. Therefore, even the results of Komiyama et al. indicate that the NW Stream can be distributed in front of M31, taking both random and systematic errors into account.

In this study, the mean metallicity of NW Stream subregions is |$[\mathrm{Fe/H}]=-1.4$|⁠. Komiyama et al. (2018b) found that the RGB and RC stars in the CMD of the NW Stream could be traced by an isochrone with |$[\mathrm{Fe/H}] = -1.37,~~~log{\mathrm{age}} = 10.00$|⁠), so the result of this study is consistent with that of previous study. Therefore, we can conclude the NW Stream is metal-poor.

Fig. 18 shows the results for the Eastern Shell and the Western Shell. The orange and red points correspond to the Eastern Shell and the Western Shell in this study. The blue, green, and yellow dots correspond to the previous studies (McConnachie et al. 2003; Ferguson et al. 2005; Tanaka et al. 2010, respectively). For both shell structures, the distances are consistent with the previous studies within the error range. Also, the metallicities of these substructures are consistent with those of the GSS. In addition to being consistent with previous studies, the results of the Western Shell and the Eastern Shell are in good agreement with each other. Simulation studies have suggested that these two shell structures were formed during the accretion of the GSS progenitor (e.g. Fardal et al. 2006, 2007). Other than these studies, this association is supported by the recent chemodynamical studies (Escala et al. 2022; Dey et al. 2023). The estimated metallicities of the two shells in this study are the almost same value, and are also consistent with those of GSS. Therefore, these results emphasize that these two shells have the same origin as GSS.

In the south-eastern part of the M31 halo, the substructures are crowded (e.g. GSS, Stream C, Stream D, Eastern Extent), so it is difficult to separate them based on photometric data only. However, there are no differences between parts of the same structure. Therefore, we think that we focus on the same substructure between adjacent subregions without the effect of the other structures, when we see within certain substructure.

3.2 The global halo properties

In this section, we use the NRGB stars defined in Section 2.5 to show the metallicity distribution in Section 3.2.1 with its estimation method, and the surface brightness profile in Section 3.2.2. Note that NRGB stars are cut by the colour–magnitude information in addition to |$p_{\mathrm{M31}} \gt 0.9$|⁠. In particular, it should be noted that the results shown in Section 3.2.1 are biased against metal-rich stars because stars with |$({\it g}-{\it i})_0\lt 2.5$| are selected.

In some previous studies (e.g. Gilbert et al. 2012; Ibata et al. 2014), metallicity distribution and surface brightness profiles were constructed by masking the substructures. However, it is difficult to mask the substructure with only the data used in this study. Therefore, it should be noted that the global halo profiles include substructures.

3.2.1 Metallicity distribution

We apply the method described in Section 2.7 to estimate the photometric metallicity of the individual stars. As described in Section 2.5, there are remaining dwarf stars with high |$p_{\mathrm{M31}}$|⁠. Assuming the distance of M31, these stars are estimated to be metal-rich stars on the CMD. Therefore, we use NRGB stars to investigate the metallicity distribution of the M31 stellar halo without uncertainties due to foreground contaminations.

The estimated metallicity distribution of NRGB stars in the entire region is shown in Fig. 19. In Fig. 19, the black line is the cumulative distribution, the solid red line is the mean metallicity, and the dashed red line is the median value. The mean, median, and standard deviation of the metallicity for the entire analysis area are estimated to be |$[\mathrm{Fe/H}]_{\mathrm{mean}} = -0.83$|⁠, |$[\mathrm{Fe/H}]_{\mathrm{median}} = -0.72$|⁠, and |$\sigma _{[\mathrm{Fe/H}]} = +0.54$|⁠, respectively. In this study, we truncate red stars with |$({\it g}-{\it i})_0\gt 2.5$|⁠, when selecting the RGB-box (see Section 2.5) to create the metallicity distribution. Therefore, it should be noted that we may have underestimated the metal-rich side of our metallicity distribution.

Fig. 20 shows the metallicity distribution of NRGB stars in 4 annuli (projected radius |$\mathrm{r_{proj}} = 0 - 2$| deg, |$2-4$| deg, |$4-6$| deg, and |$6-10$| deg). The photometric metallicity results for each annulus are summarized in Table 4. Table 4 shows that the outer annulus is more metal-poor, confirming the spatial metallicity gradient in the M31 stellar halo (Kalirai et al. 2006). In the observed fields, the outer region is closer to the Galactic plane, which is more susceptible to contamination of the Galactic disc stars because it is known that the number of such stars increases with an exponential function (Tanaka et al. 2010; Martin et al. 2013). Galactic disc stars are generally identified as red stars (metal-rich populations assuming a distance of M31), so the metallicity distribution and metallicity profile of the M31 halo would have been estimated as metal-rich, especially in the north-western region, unless the contamination sources are correctly removed. However, as shown above, the metallicity distribution in this study shows that the outer region of the M31 halo is more metal-poor than the inner region. This result indicates that NB515 is a powerful tool to eliminate foreground stars and that we are able to capture the true structure of the M31 stellar halo.

To further obtain the global metallicity property of the M31 stellar halo, we construct both the radial metallicity profile and the profiles along the south-east and north-west directions of the minor axis of the M31 disc (hereafter, we call these two profiles SE/NW minor axis profiles). In general, the minor axis profiles allow us to obtain the stellar halo profile with minimal influence of both the disc component and the stellar populations in a particular direction than the radial profile, because our analysis region is anisotropic. We construct three profiles using the NRGB stars as follows. First, the M31 halo is divided into small regions. In the case of the radial profile, the stars are divided into 3|$^\prime$| areas from the M31 centre throughout the analysis region: in the case of the minor axis profiles, the SE minor axis profile is constructed by the data in the south-eastern region (PFS_FIELD_19—PFS_FIELD_33; see Table 1) and the NW minor axis profile is from the data north-western region (PFS_FIELD_1—PFS_FIELD_16 and M31_003, 004, 009, 022, 023; see Table 1). Using these data, the observed fields are divided into |$3^\prime \times 1^{\circ }$| areas, to construct the minor axis metallicity profiles. Finally, we calculate the mean/median metallicity of the sources in each divided region. Then, the error of the mean/median metallicity in each region is the standard deviation of the metallicity of the objects. It should be the variation in the stellar metallicity within a region divided by the number of stars, but this value is sufficiently small because there are a sufficient number of sources within each region.

Fig. 21 shows the radial metallicity profile of the M31 halo and Fig. 22 is the SE and NW minor axis profiles. In each panel, the black dots are the mean/median metallicity values in each binned region. We performed MCMC linear regression on these data using emcee with |$chains=60, step=110000$|⁠. The solid red lines in Figs 21 and 22 are 1000 linear regression lines randomly sampled from the posterior distribution. The slope of the regression line for the radial profile for mean (median) metallicity is |$-0.004^{+0.001}_{-0.001}$| dex kpc|$^{-1}$| (⁠|$-0.0044^{+0.0004}_{-0.0004}$| dex kpc|$^{-1}$|⁠), indicating that the halo of M31 has a metallicity gradient.

The SE minor axis profile for has a slope in the mean/median metallicity of |$0.002^{+0.001}_{-0.001}$| dex kpc|$^{-1}$| (the reason why the slope value is a positive value is that the distance from the centre of M31 is set to a negative value), and we can see the metallicity gradient in this profile. Moreover, the NW minor axis profile has a slope in the mean (median) metallicity of |$-0.003^{+0.001}_{-0.001}$| dex kpc|$^{-1}$| (⁠|$-0.0032^{+0.0005}_{-0.0005}$| dex kpc|$^{-1}$|⁠). Therefore, in the present analysis, we confirm the existence of a metallicity gradient in both radial and minor axis directions. It should be noted that substructures of M31 are included in the analysis, so the results reflect both smooth and substructure components in M31. It is difficult to separate the M31 substructures from the smooth halo using only the observational data in this study. Further kinematic information from future spectroscopic observations, like PFS, will be important, so that we can explore the metallicity properties with respect to each of the components.

3.2.2 Surface brightness profile of all resolved RGB stars

We construct the radial and two minor axis surface brightness profiles using the NRGB stars defined in Section 2.5. It should be noted that the profile we build in this study is ‘the surface brightness of all resolved RGB stars’, but hereafter we call it simply as ‘the surface brightness profile’. The surface brightness profile is constructed by adding the fluxes of the individual halo stars as follows. First, the M31 halo is divided into the same small regions as the metallicity profiles (Section 3.2.1). Second, we divide the objects in each subregion by photometric metallicity (Section 3.2.1). The data sets are divided according to (Ibata et al. 2014) |$-2.5 \lt [\mathrm{Fe/H}] \lt -1.~-1.7 \lt [\mathrm{Fe/H}] \lt -1.1,~-1.1 \lt [\mathrm{Fe/H}] \lt 0.0$|⁠. Third, the fluxes of the sources within each divided region are added and their V-band magnitude is derived according to the photometric transformation of Ibata et al. (2007). Finally, the surface brightness profile is constructed by dividing the sum of the added fluxes for each region by the area of the region. The constructed radial profile is shown in Fig. 23, and minor axis profiles are shown in Fig. 24. In addition to the surface brightness profiles for each metallicity range, the profiles for all metallicity ranges (⁠|$-2.5 \lt [\mathrm{Fe/H}] \lt 0.0$|⁠) are also included in Figs 23 and 24 for reference.

In the left panel of Fig. 24, the surface brightness tends to decrease as the distance from the centre increases, while the profiles for distances |$\sim -3 ^{\circ }$| (−40 kpc) and |$\sim -4 ^{\circ }$| (−55 kpc) show a hump in the surface brightness with |$-1.7 \lt [\mathrm{Fe/H}] \lt -1.1$|⁠. This reflects the presence of Stream C and Stream D in this region. Similarly, for |$-2.5 \lt [\mathrm{Fe/H}] \lt -1.7$| and |$-1.7 \lt [\mathrm{Fe/H}] \lt -1.1$| in the right panel of Fig. 24, the surface brightness profile is enhanced due to the presence of the low-metallicity NW Stream at the radius of |$\sim 7 ^{\circ }$| (95 kpc) from the centre. The surface brightness also increases for |$-1.1 \lt [\mathrm{Fe/H}] \lt 0$| at |$\sim 7 ^{\circ }$| (95 kpc) in the right panel, which is due to a remnant of foreground stars in the Galaxy in this region.

We then perform a model fitting to each of these profiles. It has been shown that the surface brightness in the M31 halo can be expressed as a power law (⁠|$\Sigma (R) \propto R^{-\alpha }$| or an exponential law (⁠|$\Sigma (R) \propto \exp {(-R/R_s)}$|⁠, where |$R_s$| is the scale length. e.g. Guhathakurta et al. 2005; Irwin et al. 2005; Ibata et al. 2007; Tanaka et al. 2010). Figs 23 and 24 show the power-law and exponential-law models fitted to the observed data with dotted and solid lines, respectively. Table 5 shows the results of the fitting and previous studies. For both minor axis and radial profiles, we are able to fit power-law models with power indices ranging from |$\sim -1.3$| to |$\sim -3.0$|⁠. In Pritchet & van den Bergh (1994), the M31 halo is represented by a single de Vaucouleurs’s law (⁠|$\ln I \propto -R^{1/4}$|⁠) between 0.2 and 20 kpc. The results of this study show that the power index is |$\sim -1.3$| to |$\sim -3.0$|⁠, which is different from the structure of the inner halo represented by de Vaucouleurs’s law with a power index of |$1/4$|⁠. This may reflect the fact that the M31 halo is not a single structure.

Our model fitting using the exponential law results in a scale length range from |$\sim 10$| to |$\sim 40$| kpc. Table 5 shows that the scale lengths in this study are consistent with those of Tanaka et al. (2010). However, in Figs 23 and 24, it can be seen that the exponential law reproduces the surface brightness near the centre, but does not reproduce the M31 halo in the outer region (⁠|$\gt 80$| kpc). Therefore, the power law is more suitable than the exponential law for reproducing the entire M31 halo for the radial range of this study.

The higher the metallicities of the samples in the profiles, the smaller the |$\alpha$| in the power law and the smaller the scale length in the exponential law, indicating a steeper distribution for the metal-rich objects. This result suggests the fact that more massive (metal-rich) mergers sink further into the halo and the outer halo consists of mostly minor (metal-poor) mergers. In addition, this also is consistent with hierarchical |$\Lambda$|CDM paradigm where the stellar halo is predicted to have been mainly formed by accretion events (e.g. Johnston et al. 2008; Sharma et al. 2011). In these studies, the metal-rich haloes are formed by massive accretion events. More massive the accreting system, the stronger the dynamical friction, so massive accretion events dominate the inner halo and will have a sharper brightness profile. In contrast, the outer halo is dominated by a lot of small accretion events and should have a shallower brightness profile. Instead of this scenario, the origin of these metallicity gradients might be the gradients in the accreted dwarf galaxies. This result is also confirmed by the metallicity estimation results in Section 3.2.1.

The power indices of the radial and minor axis profiles with low metallicity (e.g. |$-2.5\lt [\mathrm{Fe/H}]\lt -1.7$|⁠) in this study do not agree with those of Ibata et al. (2014). In particular, the profiles with |$-2.5\lt [\mathrm{Fe/H}]\lt -1.7$| in this study tend to be shallower than that in Ibata et al. This may be due to the inclusion of substructures (e.g. NW Stream) in the present profiles, resulting in higher surface brightness in the outer region and shallower gradients in the profiles. In fact, Ibata et al. (2014) constructed the surface brightness profile by masking the substructure after statistically removing the contamination sources beforehand, suggesting that masking the substructure in the same way will result in a steeper gradient in our surface brightness profiles.

4 DISCUSSION

4.1 Accuracy of the M31 RGB stars selection

In Section 2.5, we calculate the probability |$p_{\mathrm{M31}}$| that the stars are likely to be RGB stars in M31 based on NB515 and Galactic latitude information. In this subsection, we discuss its accuracy of |$p_{\mathrm{M31}}$| by comparing with available spectroscopic information of the stars.

Spectroscopic observations were conducted targeting the inner halo of the M31 (e.g. Ibata et al. 2004; Chapman et al. 2008; Collins et al. 2009). Besides these pencil beam observations, some large-scale spectroscopic surveys were also conducted recently (e.g. Gilbert et al. 2009; Escala et al. 2022; Dey et al. 2023). To compare our results, we use the SPLASH data (Gilbert et al. 2018), which have many spectroscopically identified stars and are labelled as whether each star is a giant or dwarf star. Based on the line-of-sight velocity, the strength of the Na i doublet absorption line, colour, magnitude, and metallicity, the SPLASH data are attached the labels to individual stars that are likely to be a dwarf or RGB star (M31 CLASS). In this label, the positive value (⁠|$+1$|⁠, |$+2$|⁠) means that the object is likely to be an M31 RGB star, and the negative value (⁠|$-1$|⁠, |$-2$|⁠) means that it is another source such as the foreground dwarf star. Therefore, we compare this label with |$p_{\mathrm{M31}}$|⁠. Fig. 25 shows a two-dimensional histogram with M31 RGB probability on the vertical axis, and SPLASH label on the horizontal axis. In this figure, the plots in the upper right and lower left are those considered to be M31 members/non-members in both the HSC and SPLASH catalogues. To obtain the accuracy rate of our RGB selection, we evaluate it from the following fraction: the number of objects determined to be RGB/MS in HSC and SPLASH divided by the number of objects cross-matched in HSC and SPLASH.

where the RGB (dwarf) stars in the HSC catalogue are labelled as the stars with |$p_{\mathrm{M31}} \gt 0.5$| (⁠|$p_{\mathrm{M31}} \lt 0.5$|⁠), and the RGB (dwarf) stars in the SPLASH catalogue are labelled as the stars with |$\mathrm{M31CLASS} \gt 0$| (⁠|$\mathrm{M31CLASS} \lt 0$|⁠). Based on this formula, we calculate the accuracy rate of our RGB selection to be 90.2 per cent, thereby demonstrating that NB515 is a powerful tool for separating RGB and dwarf stars.

4.2 Comparison with simulation

Simulation studies for reproducing these M31 substructures have been conducted by several researchers (e.g. Fardal et al. 2007; Mori & Rich 2008; Hammer et al. 2018). In this subsection, we compare our analysis for GSS and NW Stream with some simulation results.

Many simulations have been performed to reproduce the orbits of the GSS. In particular, two hypotheses for the formation scenario of the GSS have been discussed theoretically: the minor merger hypothesis (Fardal et al. 2007; Mori & Rich 2008; Kirihara et al. 2017b), which postulates that small galaxies collided about 0.6 Gyr ago, and the major merger hypothesis (Hammer et al. 2018), which postulates that massive galaxies collided with each other a few Gyr ago. Since the distance values estimated in this study are statistically consistent with the previous studies (McConnachie et al. 2003; Conn et al. 2016), it is difficult to constrain the two hypotheses using the line-of-sight distance. However, two of the three substructures (Metal-Poor Cloud and SE Stream) detected in this study might be used to constrain the two hypotheses. For the SE Stream, similar structures have been identified in Miki, Mori & Rich (2016) and Hammer et al. (2018), so this structure is expected to form regardless of the mass of the progenitor. As for the Metal-Poor Cloud, its existence was predicted in Kirihara et al. (2017b). The details of the origin and nature of these substructures are discussed in Section 4.3, but it is suggested that two of the three structures detected in this study may need to be minor merger in order to reproduce simultaneously in the GSS formation event. However, the possibility that these structures are independent of the GSS cannot be rejected, and further kinetic and scientific information is needed.

For the NW Stream, Kirihara et al. (2017a) proposed two groups of orbits (Case A and B) that reproduce the NW Stream (see their fig. 4). In Case A, the progenitor accretes from the north-west of the M31 halo towards the centre, resulting in the distance for most parts of the NW Stream being distributed farther than that of M31. In this case, the stream shows a distance gradient where the stream is distributed far in the line-of-sight direction where |$\eta$| is large and closer where |$\eta$| is small. In contrast, in Case B, the progenitor accretes from the eastern part of the M31 halo towards the centre, and the distance of the stream is distributed entirely in front of M31. Our analysis here confirms that the NW Stream is distributed in the front of M31 where |$\eta$| is small, and in the back where |$\eta$| is large. Therefore, the trend of the distance gradient which is estimated by this study is consistent with that in Case A, although the slope of the distance gradient is slightly different. Also, in Case B, its distance distribution shows that the stream locates at the front of M31. On the other hand, in Case A, most part of the stream locates behind M31, and the stream near the centre of M31 locates at the front of M31. Therefore, the result of distance estimation in this study supports the orbital model in Case A. In any case, the photometric distance results in this study and the line-of-sight velocity data obtained by Subaru/PFS in future will allow us to place limits on stream formation events.

4.3 New substructures

In this paper, we have reported on three new substructures (Metal-Rich Cloud, Metal-Poor Cloud, and SE Stream) in the M31 halo (see Section 3.1.1). In this section, we discuss the origin of these substructures, accompanied by the results of the distance estimation.

Fig. 26 shows the distance distribution for the Metal-Poor Cloud, suggesting that this structure is distributed in front of M31. Previous PAndAS studies identified the stars that appear to be composed of the Metal-Poor Cloud (e.g. fig. 9 in Ibata et al. 2014), but these stars are treated as a smooth halo in these studies. If this structure is the component of the fully dynamically relaxed smooth halo, the results of the distance estimation should have values around the M31 distance (⁠|$\sim 785$| kpc), given that the smooth halo is a spherically symmetric structure. However, this structure is located in front of M31, and the stellar populations are concentrated in metal-poor ranges, so we conclude that this structure is not a dynamically relaxed structure and we consider that this cloud is one of the new substructures in the M31 stellar halo.

Since this region is adjacent to the Western Shell and the NW Stream, there are three possible origins for this substructure: (1) an origin from the accretion event with the GSS, (2) the NW Stream, or (3) an independent accretion origin. As is evident from Fig. 26, the photometric metallicity of this substructure is in almost agreement with the NW Stream, and also the distance distribution appears continuous from the NW Stream. However, if we look at the shape of the spatial distribution of the Metal-Poor Cloud, this overdensity shows a broad structure, unlike the thin shape of the NW Stream. Instead, the Metal-Poor Cloud has a smooth shape, which appears to be continuous from the Western Shell. Other than the spatial distribution, the standard deviations of photometric metallicities of this overdensity are larger than those of the NW Stream. Moreover, the results in Section 3.1.2 also suggest that the GSS does not only have a metal-rich population but also contains a metal-poor component, so the metallicity distribution of the GSS includes the metallicity range of the Metal-Poor Cloud. In addition to the observational evidence, according to the simulation study (Kirihara et al. 2017b), a metal-poor substructure is predicted to be outside the Western Shell originating from the GSS accretion event (called ‘Outer Western Shell’ by Kirihara et al. 2017b). Therefore, we think that it is highly likely that this structure was formed by the accretion event of the GSS. We note that it is yet uncertain from only current photometric data whether this substructure is indeed related to the formation of the GSS, the NW Stream, or individual accretion. To resolve this question, a much deeper and wider spectroscopic observation (e.g. Subaru/PFS) to measure the stellar kinematics will be needed to understand the origin of this halo structure in detail.

In addition to the Metal-Poor Cloud, we have also found the Metal-Rich Cloud at |$(\xi ,\eta) \sim (-1^{\circ },2^{\circ })$| in Fig. 11 with |$-0.6\lt [\mathrm{Fe/H}]\lt -0.1$|⁠. This structure may originate from (1) the GSS accretion event, or (2) an independent accretion event. Similar to the Metal-Poor Cloud, the metallicity range of the GSS covers the metallicity of the Metal-Rich Cloud. This overdensity is also spatially spread out and distributed near the centre of M31, where the GSS is dominant. Therefore, we think that the Metal-Rich Cloud is associated with the GSS, as is the case for the Metal-Poor Cloud.

Finally, we have found the SE Stream which is distributed in the M31 inner halo. Many of the substructures in the M31 inner halo (e.g. the Western Shell, the North-Eastern Shell, and the South-Eastern Shell) are thought to have been formed when the GSS was formed. In fig. 3 a of Miki et al. (2016) which simulated the GSS formation, stream-like structures can be seen in addition to the GSS and shell structures. Although the SE Stream is distributed outside of the locations identified in fig. 3 a of Miki et al. (2016), such a stream-like structure may have been formed during the GSS accretion, according to this simulation. Furthermore, Hammer et al. (2018) conducted the simulation for the formation of the GSS, and this simulation model (Hammer et al. 2018, see ‘#255’ model shown in fig. 8) has almost the same shaped substructure as that of the SE Stream near the centre of M31, although this simulation is a different formation scenario where the progenitor is massive. So this stream is possible to be one of the substructures associated with the GSS. Also, this stream has similar metallicity to the Stream C and D. This suggests that it may have formed further inward at the same time as the formation of the Stream C and D, and/or associated with these streams. In any case, it is not possible to determine whether the SE Stream is associated with other substructures using only the photometry in this study. Therefore, we expect that future spectroscopic observations measuring radial velocities and detailed chemical abundances of stars will allow us to elucidate the properties of this substructure, and its origin will be revealed by comparing it with more detailed simulation studies.

5 CONCLUSIONS

We have carried out a wide-area narrow-band imaging survey of the stellar halo in M31 using HSC on the Subaru Telescope. This survey covers about 50 deg|$^2$| field consisting of 33 HSC pointings in the narrow-band filter, NB515. Combining the PAndAS/g and i bands, we have constructed a method to statistically extract M31 RGB stars (NRGB stars) using the colour–colour diagram and their galactic latitude. Using this method, we have succeeded in deriving the most likely inner halo structure of M31 without suffering from as large a degree of contamination from foreground stellar populations.

The spatial distribution of NRGB stars shows prominent stream features, including already-known streams (e.g. GSS, Western Shell, NW Stream). We have also performed the distance and metallicity estimates for these substructures, and found that the properties of the GSS identified in this study are in good agreement with previous studies. Moreover, the NW Stream shows a steeper distance gradient than that derived in the previous study, suggesting that it was formed in an orbit closer to the MW. For the Streams C and D, we have resolved the distance gradients of these streams thanks to the NB515-selection and the distance and metallicity of these streams are consistent with each other. In addition to these existing structures, we have discovered three new substructures by taking advantage of NB515. Estimating photometric metallicities and metallicity dispersions, all substructures can be the tidal debris of the accretion event, which constructed the GSS.

Finally, we have derived the global halo photometric metallicity distribution and surface brightness profile using NB515-selected halo stars. The surface brightness of the metal-poor stellar population can be fitted to a shallower power-law profile with |$\alpha = -1.65 \pm 0.02$| than that of the metal-rich population with |$\alpha = -2.82 \pm 0.01$|⁠. In contrast to the relative smoothness of the halo density profile, its metallicity distribution appears to be spatially non-uniform with non-monotonic trends with radius, suggesting that the halo population had insufficient time to dynamically homogenize the accreted populations.

The accuracy of our NB515-based selection reaches 90 per cent, indicating that this narrow-band filter is a powerful tool for the analysis of existing outer halo data and deep photometry of the inner halo, which are currently under analysis (Ogami et al., in preparation). In addition, it is also expected to provide ideal targets for future large-scale spectroscopic observations.

ACKNOWLEDGEMENTS

We acknowledge support in part from MEXT Grant-in-Aid for Scientific Research (nos. JP18H05437 and JP21H05448 for MC; nos. JP21K13909 and JP23H04009 for KH). MGL was supported by the National Research Foundation grant funded by the Korean Government (NRF-2019R1A2C2084019). This work was partially supported by Overseas Travel Fund for Students (2023) of Astronomical Science Program, The Graduate University for Advanced Studies, SOKENDAI. Data analysis was in part carried out on the Multiwavelength Data Analysis System operated by the Astronomy Data Center (NAOJ/ADC) and the large-scale data analysis system co-operated by the ADC and Subaru Telescope, NAOJ. ENK acknowledges support from NSF CAREER grant AST-2233781. Based in part on data collected at Subaru Telescope and obtained from the SMOKA, which is operated by the NAOJ/ADC. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant no. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant no. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

DATA AVAILABILITY

The HSC raw images are publicly available at www.https://smoka.nao.ac.jp. The PAndAS images and reduced photometry catalogue are also available at https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/community/pandas/query.html.

Footnotes

REFERENCES

APPENDIX A: THE METHOD OF ARTIFICIAL STAR TEST

In this section, we describe how to calculate detection completeness using artificial stars in the hsciPipe (especially after version 6). In this study, we develop a python code injectstar.py to insert artificial stars into coadd-images and calculate the detection completeness using this code and hscPipe.

The general algorithm of the injectstar.py is almost the same as Aihara et al. (2018). The algorithm of the injectstar.py consists of two steps. To insert artificial stars, the first step in the injectstar.py is to obtain the point spread function (PSF) from each stacking image which is |$4000~\mathrm{pixels} \times 4000~\mathrm{pixels}$| region called patch in the hscPipe. Using these obtained PSF, the second step is to construct the artificial stars given the apparent magnitude and inject in each co-added patch image, every 60 pixel by 60 pixels. After these processes, we apply the hscPipe detection and photometry algorithm (detectCoaddSources.py and multiBandDriver.py) to the images with embedded artificial stars. The series of these processes are repeated in 0.25 mag over a range of magnitudes |$18 \le {\it NB515} \le 28$| to obtain the detection completeness, the percentage of detected artificial stars, for all patches.

Fig. A1 shows the patch images, which are before and after performing the injectstar.py. The magnitude of the embedded artificial star is 17 mag, and it can be seen that these stars are equally embedded in this patch image.