Testing the cosmological principle with the Pantheon+ sample and the region-fitting method

The $\Lambda$CDM model is generally accepted as the standard cosmological model, which is consistent with most astronomical observations (Scolnic et al. 2018; Abbott et al. 2019; Khadka & Ratra 2020; Brout et al. 2022; Cao & Ratra 2022; Dainotti et al. 2022b; Jia et al. 2022; Liu et al. 2022a; Porredon et al. 2022; Wang et al. 2022; de Cruz Pérez et al. 2023; Khadka et al. 2023; Li et al. 2023b). It is based on the fundamental assumption of the cosmological principle, namely that the Universe is statistically isotropic and homogeneous on sufficiently large scales. Despite its many successes, there have been also several analyses of observations that indicate that the Universe may be inhomogeneous and anisotropic; for instance, the fine-structure constant (Webb et al. 2011; King et al. 2012; Li & Lin 2017; Milaković et al. 2022), the direct measurement of the Hubble parameter (Bonvin et al. 2006; Koksbang 2021), the cosmic microwave background (CMB) (Bennett et al. 2003; Tegmark et al. 2003; Bielewicz et al. 2004; Bennett et al. 2011; Gruppuso et al. 2011; Ghosh & Jain 2020; Planck Collaboration et al. 2020c), the anisotropic dark energy (Koivisto & Mota 2008; Mariano & Perivolaropoulos 2012; Bayron Orjuela-Quintana et al. 2020; Motoa-Manzano et al. 2021), the large dipole of radio source counts (Rubart et al. 2014; Colin et al. 2017; Rameez et al. 2018; Singal 2019, 2023), quasar dipoles (Hutsemékers et al. 2005; Pelgrims & Hutsemékers 2016; Tiwari & Jain 2019; Hu et al. 2020; Secrest et al. 2021; Zhao & Xia 2021a, b; Dam et al. 2023), the anisotropic Hubble constant (Luongo et al. 2022; McConville & Colgáin 2023) and the type Ia supernovae (SNe Ia) dipole (Colin et al. 2011; Yang et al. 2014; Javanmardi et al. 2015; Sun & Wang 2018; Wang & Wang 2018; Tang et al. 2023), and local matter underdensity (Kazantzidis & Perivolaropoulos 2020). These analyses hint that the Universe may have a local void and a preferred expanding direction. In addition, it is assumed that the $\Lambda$CDM model also triggers a serious Hubble constant discrepancy between the Planck CMB (Planck Collaboration et al. 2020b) and the local distance ladder (Riess et al. 2019, 2022). This is known as the Hubble tension, and its statistical significance has reached 5.0 $\sigma$. Such a high confidence level could not be explained by systematic uncertainty alone, and might imply new physics beyond the $\Lambda$CDM model. There has been intense discussion focused on this issue, and a lot of theoretical explanations have been proposed. We refer the readers to some review articles (Di Valentino et al. 2021; Shah et al. 2021; Abdalla et al. 2022; Perivolaropoulos & Skara 2022; Hu & Wang 2023; Kumar Aluri et al. 2023; Kroupa et al. 2023; Riess & Breuval 2023; Vagnozzi 2023) for more detailed information about the cosmological anomalies and tension.

It is worth noting that some researchers claim that considering a void model can successfully explain the cosmic dipole and the Hubble tension (Böhringer et al. 2020; Haslbauer et al. 2020; Luković et al. 2020; Camarena et al. 2022; Cai et al. 2022; Mohammadi et al. 2023). Haslbauer et al. (2020) showed that the KBC void (Keenan et al. 2013) could naturally resolve the Hubble tension in Milgromian dynamics for the first time. By computing the Hubble constant in an inhomogeneous universe and adopting model selection via both the Bayes factor and the Akaike information criterion, Camarena et al. (2022) found that the lambda Lemaître-Tolman-Bondi ($\Lambda$LTB) model is favored with respect to the $\Lambda$CDM model at low-redshift range (0.023 $<z<$ 0.15), and this can be used to explain the Hubble tension. After that, Cai et al. (2022) proposed that a gigaparsec-scale void can reconcile the CMB and quasar dipolar tension. If considering a large and thick void, their setup can also ease the Hubble tension. At the same time, there are some findings that could be explained by a local void model. For example, inspired by the H0LiCOW results (Millon et al. 2020; Wong et al. 2020), Hu & Wang (2022b) reported a late-time transition of $H_{0;}$ that is, $H_{0}$ changes from being consistent with the CMB result to being consistent with the distance ladder one from an early-to-late cosmic time, which can be explicated by the local void. The late-time transition of $H_{0}$ was found from the observational Hubble parameter $H(z)$ data combining the Gaussian process (GP) method (Pedregosa et al. 2011) and can be used to effectively relieve the Hubble tension (a mitigation level of around 70 %). In addition, a similar $H_{0}$ descending behavior has also been discovered by utilizing various observations (such as SNe Ia, $H(z)$, baryon acoustic oscillations and megamasers) and their combinations (Krishnan et al. 2020, 2021b; Dainotti et al. 2022a; Ó Colgáin et al. 2022; Horstmann et al. 2022; Colgáin et al. 2022; Jia et al. 2023; Malekjani et al. 2023). Of course, there are also some opposing voices believing that the void model alone cannot solve the Hubble tension (Kenworthy et al. 2019; Cai et al. 2021; Castello et al. 2022). Therefore, further research on this controversial topic is necessary and would certainly be worthwhile. So far, there has been no research using the Pantheon+ sample to simultaneously map matter-density ($\Omega_{m}$) distribution and the Hubble expansion ($H_{0}$) distribution to test the cosmological principle.

In this work, we tested the cosmological principle by the region fitting (RF) method with the latest Pantheon+ sample. Compared to the previous work (Krishnan et al. 2021a, 2022; McConville & Colgáin 2023), we improved our research methodology and carried out the necessary statistical analyses. Usually, considering the simple flat $\Lambda$CDM model, $\Omega_{m}$ and $H_{0}$ are considered to be negatively correlated. Therefore, for the convenience of analysis, one parameter is usually fixed. For example, some recent works fix $\Omega_{m}$ to 0.30 and regard $H_{0}$ as a free parameter (Krishnan et al. 2021a, 2022; McConville & Colgáin 2023). In our fitting calculation process, all cosmological parameters ($\Omega_{m}$ and $H_{0}$) are free to investigate the local properties of our Universe. We would like to plot the all-sky distributions of $\Omega_{m}$ and $H_{0}$ to find out the local matter underdensity region and the preferred direction of expansion (cosmic anisotropy), respectively. The influence of redshift and the screening angle $\theta_{\mathrm{max}}$ on the final results is also considered. We found the suitable angle of the RF method for the Pantheon+ sample. Then, we analysed a combination of the local void (Enqvist 2008; Garcia-Bellido & Haugbølle 2008; Keenan et al. 2013; Wang & Dai 2013; Hwang et al. 2016; Hamaus et al. 2022; Shim et al. 2023), cosmic anisotropy (Sun & Wang 2019; Akarsu et al. 2020; Migkas et al. 2020, 2021; Akarsu et al. 2022; Horstmann et al. 2022; Rahman et al. 2022; Akarsu et al. 2023; Dhawan et al. 2023; Ebrahimian et al. 2023), and Hubble tension (Chen 2019; Riess et al. 2019; Verde et al. 2019; Planck Collaboration et al. 2020b; Riess 2020; Riess et al. 2022; Capozziello et al. 2023) in detail. Finally, we compared our results with those of previous similar research (Antoniou & Perivolaropoulos 2010; Cai & Tuo 2012; Kalus et al. 2013; Wang & Wang 2014; Yang et al. 2014; Chang & Lin 2015; Lin et al. 2016b; Chang & Zhou 2019; Hu et al. 2020; Luongo et al. 2022; McConville & Colgáin 2023) and other observations, including the CMB dipole (Planck Collaboration et al. 2016, 2020a), dark flow (Abdalla et al. 2022), bulk flow (Turnbull et al. 2012; Feix et al. 2017; Watkins et al. 2023), and galaxy cluster (Migkas et al. 2021).

The outline of this paper is as follows. In Sect. 2, we give a detailed description of the Pantheon+ sample, including redshift distribution, location distribution, and corresponding density contour, and compare it to the Pantheon sample. Section 3 briefly introduces the RF method, which we used to map the all-sky distribution of cosmological parameters. In Sect. 4, we present the results from the whole and low-redshift Pantheon+ sample ($z<0.30$, lp+) and discuss the impact of RF method with different screening angles $\theta_{\mathrm{max}}$ on the final results from the Pantheon+ sample. The corresponding investigation and discussion are given in Sect. 5. Finally, conclusions and perspectives are presented in Sect. 6.

Pantheon+, as the latest sample of SNe Ia, consists of 1701 SNe Ia light curves observed from 1550 distinct SNe and covers redshift range from 0.001 to 2.26 (Brout et al. 2022; Scolnic et al. 2022). The redshift distributions of the Pantheon (Scolnic et al. 2018) and Pantheon+ samples are shown in Fig. 1. From the redshift distribution of these two samples, it is not difficult to find that there are two main differences between the new sample and the Pantheon sample. One is that the SN number has increased significantly at low redshift. There are more than 700 additional SNe in the range of $z<$ 0.8, of which more than 500 SNe originated from $z<$ 0.08. This is mainly because the new sample adds five large samples, including the Foundation Supernova Survey (Foundation; Foley et al. 2018), the Swift Optical/Ultraviolet Supernova Archive (SOUSA; Brown et al. 2014), the Lick Observatory Supernova Search (LOSS1; Ganeshalingam et al. 2010), the second sample from LOSS (LOSS2; Stahl et al. 2019), and the Dark Energy Survey Year 3 (DES; Brout et al. 2019; Smith et al. 2020), all of which are low-redshift surveys, except DES. The other difference is that there is a significant deletion in Pantheon+ statistics between $0.8<z<1.0$. The reason is that Scolnic et al. (2022) did not use SNe from the Supernova Legacy Survey (SNLS) at $z>$ 0.8 considering sensitivity to the $U$ band in model training (56 SNe in total).

Figure 2 shows the distribution (left panel) and corresponding density (right panel) of the Pantheon+ sample in the sky of the galactic coordinate system. Some previous work has pointed out that the Pantheon SNe Ia are not uniformly distributed in the sky; half of them are located in the galactic southeast. Some SNe are very concentrated, forming a belt-like structure that is the SDSS sample (Zhao et al. 2019). At the same time, Zhao et al. (2019) also investigated the effect of the inhomogeneous distribution of the Pantheon sample on the cosmic anisotropy and found that the belt-like structure plays the most important role in the Pantheon sample. In the left panel, to highlight the difference between these two samples, we mark the newly added SNe Ia in red. From the distribution of the Pantheon+ sample, we find that the distribution of the new data is still uneven across the sky. There are only fewer observations near (l, b)$\sim$(300^∘, -30^∘). Its incomplete sky coverage is primarily due to the fact that the Pantheon+ sample consists of different subsamples, following multiple measurement strategies. In order to make it easier to comprehend this focus, we also plot the density distribution of the Pantheon+ sample utilizing the $plt.contour()$ function¹¹1https://matplotlib.org/stable/gallery/images_-contours_-and_-fields/ir regulardatagrid.html, as shown in the right panel of Fig. 2. Color-coded values represent sample fraction per unit area. From the density distribution, we can more intuitively feel the inhomogeneity of sample distribution. As can be seen, the belt-like part of the Pantheon+ sample still plays the most important role, as in the case of the Pantheon sample, while the maximum density value changes from 0.24 to 0.21 (Hu et al. 2020). This means that the Pantheon+ sample is more uniform than the Pantheon sample. The newly added SNe Ia effectively weaken the dominance of the belt structure in the distribution. The Pantheon+ sample that is relatively uniform and rich in low redshift is very suitable for analyzing the local property of the Universe (Andrade et al. 2018; Luongo et al. 2022; Kalbouneh et al. 2023). Here, we used the RF method combined with the Pantheon+ sample to describe the all-sky distribution of cosmological parameters to study the local structure of the Universe, and the detailed analysis process is shown in the next section.

The hemisphere comparison (HC) method proposed by Schwarz & Weinhorst (2007) has been widely used in the investigation of the cosmic anisotropy, such as the anisotropy of cosmic expansion (Deng & Wei 2018a; Zhao & Xia 2022), the acceleration scale of modified Newtonian dynamics (Zhou et al. 2017; Chang et al. 2018; Chang & Zhou 2019), and the temperature anisotropy of the CMB (Hansen et al. 2004; Bennett et al. 2013; Ghosh et al. 2016; Ferreira & Quartin 2021). The RF method we used is similar to it. Here, we provide a detailed introduction to this method. Its goal is to map the all-sky distribution of the cosmological parameters. The most important step is to generate random directions $\hat{D}$ $(l$, $b)$ to pick out SNe data located in specific regions to construct subdatasets, where $l\in(0^{\circ}$, $360^{\circ})$ and $b\in(-90^{\circ}$, $90^{\circ})$ are the longitude and latitude in the galactic coordinate system, respectively. The specific region is given by condition ($\theta<\theta_{\mathrm{max}}$). $\theta$ is the angle between the random direction $\hat{D}$ $(l$, $b)$ and the SN position. Here we refer to $\theta_{\mathrm{max}}$ as the screening angle used to obtain regions of different sizes, and the value range is (0^∘, 180^∘). The remaining steps can be divided into three parts, which we present in the following subsections.

We first give the best fitting results in the flat $\Lambda$CDM model employing the full Pantheon+ sample, $\Omega_{m}$ = 0.36$\pm$0.02, and $H_{0}$ = 72.83 $\pm$ 0.23 km s^-1 Mpc^-1. The results are in line with the previous research (Brout et al. 2022), except that the 1$\sigma$ error of $H_{0}$ is reduced. The main reason is that we utilized the standardized distance modulus ($\mu_{\mathrm{obs}}$; Tripp 1998), where fiducial $M$ has been determined from SH0ES 2021 Cepheid host distances (Riess et al. 2022).

After that, using the RF method with a screening angle $\theta_{\mathrm{max}}=90^{\circ}$, we mapped the all-sky distribution of $\Omega_{m}$ and $H_{0}$ and draw the 1$\sigma$ regions of local matter underdensity and cosmic anisotropy, as shown in the upper panel and lower panel of Fig. 3, respectively. In Fig. 3, the range of $\Omega_{m}$ and $H_{0}$ are (0.29, 0.41) and (72.03, 74.26), respectively. The corresponding differences are $\Delta\Omega_{m}$ = 0.12 and $\Delta H_{0}$ = 2.23 km s^-1 Mpc^-1. The values of $D_{\mathrm{max}}$ to $\Omega_{m}$ and $H_{0}$ are $D_{\mathrm{max},\Omega_{m}}$ = 3.29$\sigma$ and $D_{\mathrm{max},H_{0}}$ = 4.48$\sigma$, respectively. In the upper panel of Fig. 3, the minimum constraint of $\Omega_{m}$ is $\Omega_{m,\mathrm{min}}$ = 0.29${}^{+0.03}_{-0.02}$ (1$\sigma$) and the corresponding constraint for $H_{0}$ is 74.11${}^{+0.40}_{-0.40}$ km s^-1 Mpc^-1. The direction and corresponding 1$\sigma$ range that can be adopted to describe the local matter underdensity are

The lower panel of Fig. 3 shows the corresponding all-sky distribution of the Hubble expansion. The maximum constraint of $H_{0}$ is $H_{0,\mathrm{max}}$ = 74.26$\pm$ 0.39 km s^-1 Mpc^-1, and the corresponding constraint of $\Omega_{m}$ is 0.30${}^{+0.03}_{-0.03}$. Its confidence contours are shown in Fig. 4, marked with blue lines. The preferred direction of cosmic anisotropy and corresponding 1$\sigma$ range are

In Fig. 4, we also give the best fitting results of opposite directions; that is, $H_{0,\mathrm{min}}$ = 72.14${}^{+0.27}_{-0.27}$ km s^-1 Mpc^-1 and $\Omega_{m}$ = 0.40${}^{+0.02}_{-0.02}$ , which are marked with red lines.

The statistical isotropic results shown in Fig. 5 can be well described by Gaussian functions. For the isotropy analysis, the statistical significances ($\alpha$) of the real data are 2.78$\sigma$ for $\Omega_{m}$ anisotropy (upper, purple panel) and 3.96$\sigma$ for $H_{0}$ anisotropy (upper, blue panel). The statistical significance ($\beta$) of the real data given by the isotropy RP analysis are 2.34$\sigma$ for $\Omega_{m}$ anisotropy (lower, purple panel) and 3.15$\sigma$ for $H_{0}$ anisotropy (lower, blue panel).

All-sky distributions of cosmological parameters ($\Omega_{m}$ and $H_{0}$) from the total sample show that there is an obvious local underdensity area and a preferred direction of cosmic expansion. From Fig. 3, we can find that there are some weird features that may result from fluctuations in the matter density, or fluctuations in the Hubble expansion (Gurzadyan, V. G. et al. 2023). In the upper panel, the red and blue areas represent areas of higher and lower material density, respectively. In the lower panel, the red areas represent regions of higher Hubble expansion. The continuity of these structures suggests that the results may not be sensitive to individual SNe Ia. Combining these two panels of Fig. 3, we find that the local underdensity in the upper panel overlaps with the higher Hubble expansion (cosmic anisotropy) in the lower panel. This may imply that local matter density might be responsible for the anisotropy of the accelerated expansion of the Universe. The corresponding 1$\sigma$ regions are quite obvious. In other words, the 1$\sigma$ range of $\Omega_{m}$ is significantly larger than that of $H_{0}$. The main reason for this phenomenon might be that these two parameters have different sensitivities on the redshift. The inhomogeneous matter-density distribution could also affect the decelaration parameter $q_{0}$, causing a larger $q_{0}$ in the local underdensity. Theoretically, the maximum $q_{0}$ direction should be consistent with the local underdensity direction. For the determination of $H_{0}$ value with local distance indicators, the observed $H_{0}$ values depend on the average matter density within the distance range covered (Böhringer et al. 2020). Combining the $\Omega_{m}$ distribution, we can find that there might be a region of low matter density, which leads to a smaller average matter density which makes the $H_{0}$ measurements higher. Therefore, the local underdensity can exacerbate the magnitude of the Hubble Tension (Böhringer et al. 2020; Hu & Wang 2022b; Jia et al. 2023; Hu & Wang 2023). In other words, accounting for the local underdensity or cosmic anisotropy could alleviate the current Hubble tension.

From Fig. 5, we find that the isotropy analysis shows an obvious statistical significance, especially the study of $H_{0}$, which is nearly 4$\sigma$. Afterwards, considering the spatial inhomogeneity of the real sample, we performed the isotropy RP analysis and find a slight decrease in statistical significance. This suggests that inhomogeneous spatial distribution of a real sample can increase the deviation from isotropy. In addition, combining all the statistical results, we find the statistical significations of $H_{0}$ anisotropy are more obvious than that of $\Omega_{m}$ anisotropy. From the distributions of a single parameter, we find an obvious anisotropy. But the effect of one on the other might be canceled out, making the all-sky distribution of luminosity distances isotropic. Therefore, in order to clear up this doubt, we plotted the all-sky distribution of luminosity distances when $z$ is set to 2.26 and give the relationships between $d_{L}$ and $z$ of two anisotropic hemispheres, as shown in Figs. 14 and 15, respectively. The results of these two figures show that the all-sky distribution of luminosity distances is also anisotropic and the $d_{L}-z$ relationships obtained from two anisotropic hemispheres have obvious differences. The direction of larger luminosity distance is consistent with the anisotropy direction from the Pantheon+ sample. All in all, our investigations from the Pantheon+ sample display an inhomogeneous and anisotropic universe, and the statistical results show no low confidence level.

Motivated by recent research (Krishnan et al. 2020, 2021b; Dainotti et al. 2021, 2022a; Colgáin et al. 2022; Hu & Wang 2022b; Ó Colgáin et al. 2022; Jia et al. 2023; Malekjani et al. 2023) hinting that $H_{0,z}$ might evolve with redshift, we conducted reanalyses using the low-redshift (lp+) sample. If $H_{0,z}$ evolves with redshift (may be caused by the local void, and so on), then the measured $H_{0}$ of low-redshift and high-redshift SNe samples are different. In that case, if the redshift distribution is not uniform across the sky, this could create biases in the anisotropy detection, especially when the full sample is used; moreover, the redshift range is wider. The value of $H_{0}$ obtained from the higher redshift hemisphere is smaller than that obtained from the lower redshift hemisphere. Severely uneven redshift distribution might lead to an increase in the degree of anisotropy or a bias in the anisotropy detection. The reanalysis results are similar to those of the Pantheon+ sample and verify our previous findings. From the results of the lp+ sample (Fig. 7), we find that directions of the local matter underdensity and cosmic anisotropy given by the lp+ sample are inconsistent. Comparing the results in Figs. 4 and 8, we can find that the best fits of the two anisotropic hemispheres obtained from the lp+ sample are significantly worse than the results of the total sample. In Fig. 16, we also show the relationship between $D_{\mathrm{max}}$ and $z_{\mathrm{max}}$, the detailed information are display in Table 2. From Fig. 16, we can find that $D_{\mathrm{max}}$ changes with $z_{\mathrm{max}}$, and has a maximum value near $z_{\mathrm{max}}$ = 0.30. This might imply that the structure of the Universe changes with redshift. In addition, we also find that the influence of redshift on the all-sky distribution of $\Omega_{m}$ is larger than that of $H_{0}$. Finally, from the statistical results of lp+ sample, we also find that the statistical significance obtained from the isotropy RP analysis are lower than that from the isotropy analysis. This finding confirms that inhomogeneous spatial distribution of a real sample can increase the deviation from isotropy.

Theoretically, the screening angle $\theta_{\mathrm{max}}$ in the RF method can be any value between 0$\degr$ and 180$\degr$. However, due to the limitations of real data in quantity, spatial part, and redshift distribution, $\theta_{\mathrm{max}}$ cannot be selected arbitrarily. In Sect. 4.2, we combine the RF method with the latest Pantheon+ sample, hoping to find a suitable screening angle $\theta_{\mathrm{max}}$ . The analysis results are shown in Fig. 19. Finally, we find that the screening angle $\theta_{\mathrm{max}}$ of the Pantheon+ sample is 60$\degr$. Then, based on the 60$\degr$ RF method, we restudied the Pantheon+ sample. The corresponding results are shown in Figs. 11, 12, and 13. From Fig. 11, we find that there is a small area (314.16$\degr$, -22.92$\degr$) with a lower $H_{0}$ value in the higher $H_{0}$ area, but this structure does not appear in the RF (90$\degr$) results. The lower $H_{0}$ area may be caused by high material density structures (e.g. dark matter halo) or statistical uncertainty. Narrowing the screening angle reduces the number of SNe constraining cosmological parameters, which may increase the uncertainty of the constraint. We performed a bootstrap resampling of the sample (ignoring the direction, repeated 2000 times) to study the dependence of the best-fit result error and $D_{\mathrm{max}}$ on the number of SNe Ia used. Here, the number of SNe ranges from 100 to 700, with a total of seven groups. From the results of Fig. 17, we can find that the average error and $D_{\mathrm{max}}$ increases as the used number of SNe decreases. The number of SNe used significantly affects the average error but does not significantly affect the $D_{\mathrm{max}}$ value. The $D_{\mathrm{max}}$ difference caused by the reduction in number is much smaller than the total $D_{\mathrm{max}}$. The preferred direction of cosmic anisotropy using the RF method with $\theta_{\mathrm{max}}$ = 60$\degr$ is in line with that using $\theta_{\mathrm{max}}$ = 90$\degr$ within a 1$\sigma$ error. However, narrowing the screening angle means that the statistical significance of the isotropy analyses are significantly reduced. Interestingly, unlike previous findings, the statistical significance of the isotropy RP analysis is actually higher than that of the isotropy one.

Overall, we find that the all-sky distributions of cosmological parameters deviate significantly from isotropy, as shown in Figs. 3, 7, and 11. All preferred directions we obtained are consistent with each other within a 1$\sigma$ range, and they are in line with previous research that traced the anisotropy of $\Omega_{m}$ and $H_{0}$ (Antoniou & Perivolaropoulos 2010; Cai & Tuo 2012; Kalus et al. 2013; Chang & Lin 2015; Hu et al. 2020) and other dipole researches (Wang & Wang 2014; Yang et al. 2014; Lin et al. 2016b; Pandey 2017; Chang & Zhou 2019; Dam et al. 2023). However, they are different from those of Luongo et al. (2022) and McConville & Colgáin (2023), which are consistent with the results of the CMB dipole (Planck Collaboration et al. 2016, 2020a). In addition, comparing with other independent observations (as shown in Table 3) including the CMB dipole (Planck Collaboration et al. 2016, 2020a), dark flow (Abdalla et al. 2022), bulk flow (Turnbull et al. 2012; Feix et al. 2017; Watkins et al. 2023), and galaxy cluster (Migkas et al. 2021), it is easy to find that the directions of the larger $H_{0}$ and the smaller $\Omega_{m}$ are not consistent with the CMB dipole (Planck Collaboration et al. 2016, 2020a), but they coincide with the bulk flow (Turnbull et al. 2012; Feix et al. 2017; Watkins et al. 2023) and the galaxy cluster (Migkas et al. 2021). To facilitate understanding, we aggregated these results with the ones we obtained, marking them on the galactic coordinate system, as shown in Fig. 18. The effect of peculiar velocities and the bulk flow on SNe Ia cosmology has already been discussed (Hui & Greene 2006; Davis et al. 2011; Betoule et al. 2014; Rameez et al. 2018). They can make a tiny shift in the best-fit cosmological parameters and the preferred direction locally (Colin et al. 2019).

In this paper, we propose the RF method for the first time and combine this method with the Pantheon+ sample to test the cosmological principle. From the matter density and the Hubble expansion distributions mapped using the RF method, we find that the all-sky distributions of cosmological parameters deviate significantly from isotropy. The corresponding distribution of luminosity distance also deviates from isotropy; that is, the $d_{L}$-$z$ relation actually diverges from region to region. Results of statistical isotropy analyses (isotropy and isotropy RP) show relatively high confidence levels: 2.78$\sigma$ (isotropy) and 2.34$\sigma$ (isotropy RP) for the local matter underdensity, 3.96$\sigma$ (isotropy) and 3.15$\sigma$ (isotropy RP) for the cosmic anisotropy. Comparing the results of statistical isotropy analyses, we find that inhomogeneous spatial distribution of real sample can increase the deviation from isotropy. The statistical significations of $H_{0}$ anisotropy are more obvious than that of $\Omega_{m}$ anisotropy. This might hint that parameter $H_{0}$ is more sensitive to cosmic anisotropy. The similar results and findings are found from reanalyses of the lp+ sample and the lower screening angle ($\theta_{\mathrm{max}}$ = 60$\degr$), but with a slight decrease in statistical significance. In addition, we find that $D_{\mathrm{max}}$ changes with $z_{\mathrm{max}}$ and has a maximum value near $z_{\mathrm{max}}$ = 0.30. The average error and $D_{\mathrm{max}}$ increase as the used number of SNe decreases. Comparing with the previous researches, we find all preferred directions we obtained are in line with that were provided by Antoniou & Perivolaropoulos (2010), Cai & Tuo (2012), Wang & Wang (2014), Yang et al. (2014), Kalus et al. (2013), Chang & Lin (2015), Lin et al. (2016b), Chang & Zhou (2019), and Hu et al. (2020), but they are not consistent with those given by Luongo et al. (2022) and McConville & Colgáin (2023).

Until now, many SNe Ia have been observed and have been widely used for cosmological applications (Hoscheit & Barger 2018; Perivolaropoulos & Skara 2021; Cowell et al. 2023; Hu & Wang 2022a; Wang 2022; Briffa et al. 2023).⁴⁴4For more recent studies on the cosmological applications employing the SNe Ia sample, see Bargiacchi et al. (2023), Cao & Ratra (2023), de Jaeger & Galbany (2023), Hu et al. (2023), Jin et al. (2023), Kumar et al. (2023), Lapi et al. (2023), Mandal et al. (2023), Ó Colgáin et al. (2023), Pastén & Cárdenas (2023), Perivolaropoulos & Skara (2023), Sakr (2023), Van Raamsdonk & Waddell (2023), Wang et al. (2023a), and references therein. However, there is still an inhomogeneous distribution in data, which significantly affects the testing of cosmological principle. One way to solve this problem is to add some new SNe Ia measurements, and another way is to consider other independent observations; for example, quasar (Lusso & Risaliti 2016; Bisogni et al. 2017; Lusso & Risaliti 2017; Risaliti & Lusso 2019; Cao et al. 2022b; Khadka et al. 2022; Liu et al. 2023a; Zajaček et al. 2023), gamma-ray burst (GRB; Wang et al. 2015; Shirokov et al. 2020; Hu et al. 2021; Cao et al. 2022a; Liu et al. 2022b; Lovyagin et al. 2022; Liang et al. 2022; Wang et al. 2022), fast radio burst (FRB; Hagstotz et al. 2022; Wu et al. 2022; James et al. 2022; Zhao et al. 2022; Gao et al. 2023; Liu et al. 2023b), Tip of the Red Giant Branch (TRGB; Freedman et al. 2019, 2020; Freedman 2021), galaxy cluster (Javanmardi & Kroupa 2017; Migkas et al. 2020, 2021), and gravitational wave (GW; Abbott et al. 2017; Chen et al. 2018; Jin et al. 2023; Wang et al. 2023b) observations, among others. At present, these observations may still have some deficiencies in quantity or quality, but this situation is expected to improve in the future. The e-ROSITA all-sky survey (Merloni et al. 2012; Predehl 2012; Kolodzig et al. 2013; Lusso 2020), Einstein Probe (EP; Yuan et al. 2015), French–Chinese satellite space-based multi-band astronomical variable objects monitor (SVOM; Wei et al. 2016), China Space Station Telescope (CSST) photometric survey (Xu et al. 2022; Miao et al. 2023; Li et al. 2023a), and Transient High-Energy Sky and Early Universe Surveyor (THESEUS; Amati et al. 2018) space missions together with ground- and space-based multi-messenger facilities will provide a lot of observations and measurements, improve the observational quality, and probe the poorly explored high-redshift universe. The Australian Square Kilometer Array Pathfinder (ASKAP; Chapman et al. 2017), MeerKAT (Sanidas et al. 2018), Very Large Array (VLA; Law et al. 2018), and Canadian Hydrogen Intensity Mapping Experiment (CHIME)/FRB Outriggers (Leung et al. 2021) will provide a large number of positioned FRBs in the future, which will give higher precision cosmological constraints. In addition, the Advanced Laser Interferometer Gravitational-wave Observatory (aLIGO; LIGO Scientific Collaboration et al. 2015) and Virgo (Acernese et al. 2015) detectors will provide more GW events. In combination with the electromagnetic counterparts, model-independent constraints will be given on the cosmological parameters. The high-quality observations enable us to examine the cosmological principle at higher redshifts and investigate whether the Hubble tension is related to the failure of the cosmological principle.

We thank the anonymous referee for constructive comments. This work was supported by the National Natural Science Foundation of China (grant No. 12273009), the China Manned Spaced Project (CMS-CSST-2021-A12), Jiangsu Funding Program for Excellent Postdoctoral Talent (20220ZB59), Project funded by China Postdoctoral Science Foundation (2022M721561), NWO, the Dutch Research Council, under Vici research programme ‘ARGO’ with project number 639.043.815, Yunnan Youth Basic Research Projects 202001AU070013 and National Natural Science Foundation of China (grant No. 12303050).