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

The cosmological principle is fundamental to the standard cosmological model. It assumes that the Universe is homogeneous and isotropic on very large scales. As the basic assumption, it must stand the test of various observations. In this work, we investigated the properties of the Pantheon+ sample, including redshift distribution and position distribution, and we give its constraint on the flat Λ CDM model: Ω m = 0.36 ± 0.02 and H 0 = 72.83 ± 0.23kms − 1 Mpc − 1 . Then, using the region fitting (RF) method, we mapped the all-sky distribution of cosmological parameters ( Ω m and H 0 ) and find that the distribution significantly deviates from isotropy. A local matter underdensity region exists toward (308 . 4 ◦ + 47 . 6 − 48 . 7 , − 18 . 2 ◦ + 21 . 1 − 28 . 8 ) as well as a preferred direction of the cosmic anisotropy (313 . 4 ◦ + 19 . 6 − 18 . 2 , − 16 . 8 ◦ + 11 . 1 − 10 . 7 ) in galactic coordinates. Similar directions may imply that local matter density might be responsible for the anisotropy of the accelerated expansion of the Universe. Results of statistical isotropy analyses including Isotropy and Isotropy with real-data positions (RP) show high confidence levels. For the local matter underdensity, the statistical significances are 2.78 σ (isotropy) and 2.34 σ (isotropy RP). For the cosmic anisotropy, the statistical significances are 3.96 σ (isotropy) and 3.15 σ (isotropy RP). The comparison of these two kinds of statistical isotropy analyses suggests that inhomogeneous spatial distribution of real sample can increase the deviation from isotropy. The similar results and findings are also found from reanalyses of the low-redshift sample (lp+) and the lower screening angle ( θ max = 60 ◦ ), but with a slight decrease in statistical significance. Overall, our results provide clear indications for a possible cosmic anisotropy. This possibility must be taken seriously. Further testing is needed to better understand this signal.

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 (ΛLTB) model is favored with respect to the Λ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(Krishnan et al. , 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 (Ω 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(Krishnan et al. , 2022;;Mc-Conville & Colgáin 2023), we improved our research methodology and carried out the necessary statistical analyses.Usually, considering the simple flat ΛCDM model, Ω 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 Ω m to 0.30 and regard H 0 as a free parameter (Krishnan et al. 2021a(Krishnan et al. , 2022;;McConville & Colgáin 2023).In our fitting calculation process, all cosmological parameters (Ω m and H 0 ) are free to investigate the local properties of our Universe.We would like to plot the all-sky distributions of Ω 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 θ 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. 2020Migkas et al. , 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 θ 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+ sample
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 Pan-theon+ 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 5HGVKLIW] 1XPEHU 3DQWKHRQ 3DQWKHRQ3OXV 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)∼(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() function1 , 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.

Region fitting method
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 D (l, b) to pick out SNe data located in specific regions to construct subdatasets, where l ∈ (0 • , 360 • ) and b ∈ (−90 • , 90 • ) are the longitude and latitude in the galactic coordinate system, respectively.The specific region is given by condition (θ < θ max ).θ is the angle between the random direction D (l, b) and the SN position.Here we refer to θ 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.

Fitting parameters
According to subdatasets obtained in the previous step, the corresponding best fits of the cosmological parameters are obtained by minimizing value of χ2 , where ∆µ is the difference between the observational distance modulus µ obs and the theoretical distance modulus µ th : For the flat ΛCDM model, the corresponding form of µ th can be written as here, z i is the peculiar-velocity-corrected CMB-frame redshift of each SN (Carr et al. 2022), m is the apparent magnitude of the source, M is the absolute magnitude, and d L is the luminosity distance expressed in megaparsec, defined in the following equation: where c is the speed of light.
The statistical (C stat ) and systematic covariance matrices (C sys ) are combined and adopted to constrain the cosmological parameters: (5) The datasets we used, µ obs, and C stat+sys, are provided by Brout et al. (2022) and can be obtained online 2 .C stat+sys includes all the covariance between SNe (and also Cepheid host covariance) due to systematic uncertainties (Brout et al. 2022).C stat mainly includes the full distance error and measurement noise.C sys can manifest in three key places in the analysis: (1) from changing aspects affecting the light-curve fitting; (2) from changing redshifts that propagate to changes in distance modului relative to a cosmological model; and (3) from changes in the simulations used for bias corrections (Brout et al. 2022).More detailed information about the covariance matrices can be found in Sect.2.2 of Brout et al. (2022).Unlike the previous analyses (Betoule et al. 2014;Colin et al. 2019), we did not introduce an intrinsic scatter.So, in the fitting process, there are only two free parameters (Ω m and H 0 ), which makes them strongly correlated.For the different subdatasets, Ω m and H 0 are fitted simultaneously.The C stat+sys we used were obtained by cropping the total covariance matrix according to the SNe Ia subsample.In this work, the minimization was performed employing a Bayesian Markov chain Monte Carlo (MCMC; Foreman-Mackey et al. 2013) method with the emcee package3 .All the fittings in this paper were obtained adopting this python package.The MCMC samples were plotted utilizing the getdist package (Lewis 2019).

All-sky distribution
During the calculation, we repeated 5000 random directions D (l, b), that is, 5000 sets of best fitting results (Ω m and H 0 ).Based on the fitting results, the all-sky distributions of Ω m and H 0 are mapped, respectively.From the distribution of Ω m and H 0 , we can obtain information concerning the local matter underdensity and cosmic anisotropy, respectively.In order to describe the degree of deviation from the cosmological principle, we define a parameter D max (σ) whose form is as follows: Here, P min and P max are the minimum value and the maximum value of the best fitting results, respectively.σ P min and σ P max are the corresponding 1σ error.The local matter underdensity direction and the preferred direction of cosmic anisotropy are marked by the corresponding location of the lowest Ω m and the largest H 0 .The corresponding 1σ regions are plotted in terms of the values of P 0, which is calculated from 1.00(σ) ⩾ P 0 − P fit where P fit , representing the lowest Ω m constraint or the largest H 0 constraint, is used to find the 1σ range of the local matter underdensity direction or the preferred direction.σ P fit represents the corresponding error.P 0 represents the constraints that are consistent with P fit within 1σ error, and σ P 0 are the corresponding 1σ values.We note that P 0 and σ P 0 are filtered from the total Ω m /H 0 fitting results depending on Eq. 7. Here, P 0 and P fit are completely independent and were obtained using different SNe subsamples.

Statistical analyses
In order to examine whether the discrepancy degree of the cosmological parameters from the Pantheon+ sample is consistent with statistical isotropy, we plan to carry out statistical isotropic analyses.To achieve this, we spread the original data set evenly across the sky.After that, we were able to obtain the D max for the isotropic dataset.Meanwhile, an additional isotropic analysis was also considered.We preserved the spatial inhomogeneity of real sample and then randomly distributed the real dataset, which randomly redistributed the distance moduli and redshift combination to real-data positions (RP) only.Given the limitations of computing time, we repeated it 500 times; this gave acceptable statistics.For convenience, we refer to these two approaches as isotropy analysis and isotropy RP analysis.

Results
We  After that, using the RF method with a screening angle θ max = 90 • , we mapped the all-sky distribution of Ω m and H 0 and draw the 1σ regions of local matter underdensity and cos-mic anisotropy, as shown in the upper panel and lower panel of Fig. 3, respectively.In Fig. 3, the range of Ω m and H 0 are (0.29, 0.41) and (72.03, 74.26), respectively.The corresponding differences are ∆Ω m = 0.12 and ∆H 0 = 2.23 km s −1 Mpc −1 .The values of D max to Ω m and H 0 are D max,Ω m = 3.29σ and D max,H 0 = 4.48σ, respectively.In the upper panel of Fig. 3, the minimum constraint of Ω m is Ω m,min = 0.29 +0.03 −0.02 (1σ) and the corresponding constraint for H 0 is 74.11 +0.40  −0.40 km s −1 Mpc −1 .The direction and corresponding 1σ 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,max = 74.26±0.39 km s −1 Mpc −1 , and the corresponding constraint of Ω 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σ range are In Fig. 4, we also give the best fitting results of opposite directions; that is, H 0,min = 72.14+0.27 −0.27 km s −1 Mpc −1 and Ω 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 (α) of the real data are 2.78σ for Ω m anisotropy (upper, purple panel) and 3.96σ for H 0 anisotropy (upper, blue panel).The statistical significance (β) of the real data given by the isotropy RP analysis are 2.34σ for Ω m anisotropy (lower, purple panel) and 3.15σ for H 0 anisotropy (lower, blue panel).and 2σ) and marginalized likelihood distributions for the parameter space (Ω m and H 0 ) in the spatially flat ΛCDM model from the SNe Ia subsamples, which corresponds to H 0,min (red line) and H 0,max (blue line).The best fitting results of preferred direction are H 0,max = 74.26+0.40  −0.39 km s −1 Mpc −1 , Ω m = 0.30 +0.03 −0.03 (blue line).The best fitting results of opposite directions are H 0,min = 72.14+0.27  −0.27 km s −1 Mpc −1 and Ω m = 0.40 +0.02  −0.02 (red line).
Fig. 5. Distribution of discrepancy degree D max in 500 simulated isotropic datasets.The upper two panels show the results of statistical isotropic analyses (isotropy).The lower two panels show the results of statistical isotropic analyses that preserve the spatial inhomogeneity of real data (isotropy RP).Purple and blue represent the statistical results of Ω m and H 0 .The black curve is the best fit to the Gaussian function.
The solid black and vertical dashed lines are commensurate with the mean and the standard deviation, respectively.The red lines represent the discrepancy degree from the real data.For the isotropy analyses, the statistical significance of the real data are 2.78σ for Ω m anisotropy and 3.96σ for H 0 anisotropy.For the isotropy RP analyses, the statistical significance of the real data are 2.34σ and 3.15σ for Ω m and H 0 anisotropy, respectively.

Reanalyses at low redshift
Recently, Hu & Wang (2022b) reported a late-time transition of H 0 ; that is, H 0 changes from a low value to a high one from early to late cosmic time by combining the GP method and H(z) data.The H 0 transition occurs at z ∼ 0.49.From other observations, a similar descending behavior of H 0 has been found adopting other methods (see Hu & Wang (2023) for a review of H 0 descending trend).Around z ∼ 0.40, H 0 starts to decrease (Jia et al. 2023).
Both the transition behavior and the descending trend of H 0 can effectively alleviate the Hubble tension.Kelly et al. (2023) gave new H 0 measurements from the gravitationally lensed SNe Ia Refsdal (Refsdal 1964).The redshifts of the lens and the source are 0.54 and 1.49 (Kelly et al. 2015), respectively.Utilizing eight cluster lens models, they inferred H 0 = 64.80+4.40 −4.30 km s −1 Mpc −1 .Then, using the two models most consistent with observations, they found H 0 = 66.60 +4.10 −3.30 km s −1 Mpc −1 .Marking the values of lensing redshift and H 0 on Fig. 3 in Hu & Wang (2022b) and Fig. 4 in Jia et al. (2023), we find that these results are in good agreement with the evolutionary behaviors of H 0 .Motivated by the H 0 special behaviors in the late-time universe, we plan to construct a low-redshift subsample based on the Pantheon+ sample to make a reanalyses.According to the previous research (Hu & Wang 2022b;Jia et al. 2023), and considering the need for a sufficient sample size, we decided to select 1218 SNe with redshift less than 0.30 in the Pantheon+ sample to construct sub-sample (named the lp+ sample).Detailed steps and results are as follows.
At first, we give the cosmological constraint in the flat ΛCDM model using the lp+ sample; that is, Ω m = 0.44±0.04 and H 0 = 72.43±0.30km s −1 Mpc −1 , which is consistent with the re-  sults from the full Pantheon+ sample and from Fig. 16 in Brout et al. (2022).For ease of comparison, we give the constraint results of the Pantheon+ and lp+ samples in The constraints corresponding to preferred directions of the local matter underdensity and cosmic anisotropy are Ω m,min = 0.25±0.05and H 0 = 73.76±0.41km s −1 Mpc −1 , and H 0,max = 74.11+0.47  −0.48 km s −1 Mpc −1 and Ω m = 0.27±0.06.In Fig. 8, we show the confidence contours in the preferred and opposite directions of the local matter underdensity.The isotropic statistical results from the lp+ sample are shown in Fig. 9.For the isotropy analysis, the statistical significance of the real data are 3.68σ for Ω m anisotropy (upper, purple panel) and 3.01σ for H 0 anisotropy (upper, blue panel).The isotropy RP analysis gives the statistical significance of the real data, which is 2.09σ for Ω m anisotropy (lower, purple panel) and 1.35σ for H 0 anisotropy (lower, blue panel).

Different screening angles θ max
Theoretically, if there are enough SNe data distributed at different redshifts in a certain direction, the constraints of the cosmological parameters in this direction can be obtained (Lin et al. 2016a;Deng & Wei 2018b;Kumar Aluri et al. 2023).Due to the limitation of the number of observational data, it is currently impossible to directly give the result of the limitation of cosmological parameters in a certain direction.Therefore, the 90 • RF method is usually used to derive the cosmological constraint of the single direction, that is the HC method (Schwarz & Weinhorst 2007).In this subsection, we try to find out the limit value  H 0,max = 74.81+0.71  −0.71 km s −1 Mpc −1 and Ω m = 0.28 +0.05 −0.04 .The statistical isotropy results show that the statistical confidences of the real data are 1.70σ for the isotropy analysis and 2.27σ for the isotropy RP analysis.
All results in this section are summarized in Table 1.Additionally, we provide the reduced chi-square χ 2 r for the anisotropic direction.

Discussion
All-sky distributions of cosmological parameters (Ω 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σ regions are quite obvious.In other words, the 1σ range of Ω 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 Ω 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.
Fig. 13.Distribution of discrepancy degree D max obtained from H 0 distribution in 500 simulated isotropic datasets.The left panel and right panel show the results of isotropy and isotropy RP analyses, respectively.The black curve is the best fit to the Gaussian function.The solid black and vertical dashed lines are commensurate with the mean and the standard deviation, respectively.The red lines represent the discrepancy degree from the real data.The statistical significances of the real data are 1.70σ for the isotropy analysis and 2.27σ for the isotropy RP analysis.
From Fig. 5, we find that the isotropy analysis shows an obvious statistical significance, especially the study of H 0 , which is nearly 4σ.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 Ω 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 investiga-tions 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(Krishnan et al. , 2021b;;Dainotti et al. 2021Dainotti et al. , 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 max and z max , the detailed information are display in Table 2. From Fig. 16, we can find that D max changes with z max , and has a maximum value near z 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- * D max represents the degree of deviation from the cosmological principle, and the larger the value, the higher the degree of deviation.
sky distribution of Ω 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 θ max in the RF method can be any value between 0 • and 180 • .However, due to the limitations of real data in quantity, spatial part, and redshift distribution, θ 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 θ max .The analysis results are shown in Fig. A.1.Finally, we find that the screening angle θ max of the Pantheon+ sample is 60 • .Then, based on the 60 • 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• , -22.92 • ) with a lower H 0 value in the higher H 0 area, but this structure does not appear in the RF (90 • ) 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 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 max increases as the used number of SNe decreases.The number of SNe used signif-icantly affects the average error but does not significantly affect the D max value.The D max difference caused by the reduction in number is much smaller than the total D max .The preferred direction of cosmic anisotropy using the RF method with θ max = 60 • is in line with that using θ max = 90 • within a 1σ 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σ range, and they are in line with previous research that traced the anisotropy of Ω 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 Mc-Conville & 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. 2016Collaboration et al. , 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 Ω 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).

Conclusions and perspectives
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σ (isotropy) and 2.34σ (isotropy RP) for the local matter underdensity, 3.96σ (isotropy) and 3.15σ (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 Ω 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 (θ max = 60  (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).4However, 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. 2019Freedman et al. , 2020;;Freedman 2021), galaxy cluster (Javanmardi & Kroupa 2017;Migkas et al. 2020Migkas et al. , 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) ) with 1σ range in other independent observations.We mark the directions given by Ω m with plus signs at different screening angles, and the directions given by H 0 are shown with triangles at different screening angles.The color represents the important result obtained in this work; black shows the result given by other independent 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).In combination with the electromagnetic counterparts, modelindependent 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.Notes: As this work was about to be completed, Perivolaropoulos (2023) used the HC method to test the cosmic isotropy of the SNe Ia absolute magnitudes from the Pantheon+ and SH0ES samples in various redshift/distance bins.They found that sharp changes of the level of anisotropy occuring at distances under 40 Mpc in the real samples.If there are enough local observations in the future, more local information about our Universe could be obtained by combining our method with the idea of Perivolaropoulos (2023).This will be pursued in future work.

Fig. 2 .
Fig. 2. Distribution and corresponding density contour in the galactic coordinate system.Left panel shows the SNe distribution of the Pantheon sample.Red and blue points represent the new added SNe and the SNe where the Pantheon sample and the Pan sample overlap, respectively.Right panel shows the corresponding density contour.
first give the best fitting results in the flat ΛCDM model employing the full Pantheon+ sample, Ω m = 0.36±0.02,and H 0 = 72.83± 0.23 km s −1 Mpc −1 .The results are in line with the previous research(Brout et al. 2022), except that the 1σ error of H 0 is reduced.The main reason is that we utilized the standardized distance modulus (µ obs;Tripp 1998), where fiducial M has been determined from SH0ES 2021 Cepheid host distances(Riess et al. 2022).

Fig. 3 .
Fig. 3. All-sky distribution of cosmological parameters utilizing the Pantheon+ sample combined with the 90 • RF method.The upper and the lower panels show the results of Ω m and H 0 , respectively.The corresponding values of D max (σ) are 3.29σ and 4.48σ, respectively.The star marks the directions of the lowest Ω m (upper panel) and the largest H 0 (lower panel), and the red circle outlines the corresponding 1σ areas.The directions and 1σ areas are parameterized as Ω m,min (308.4 •+47.6 −48.7 , −18.2 •+21.1 −28.8 ) and H 0,max (313.4 •+19.6 −18.2 , −16.8 •+11.1 −10.7 ).

Fig. 7 .
Fig. 7. All-sky distribution of cosmological parameters utilizing the z < 0.30 part of Pantheon+ sample (lp+ sample) combined with the 90 • RF method.The upper and the lower panel show the results of Ω m and H 0 , respectively.The corresponding values of D max (σ) are 4.17σ and 4.49σ, respectively.The star marks the directions of the lowest Ω m (upper panel) and the largest H 0 (lower panel), and the red circle outlines the corresponding 1σ areas.The directions and 1σ areas are parameterized as Ω m (26.4 •+26.3 −100.9 , −19.3 •+35.4 −17.0),H 0 (321.9•+72.5 −33.5 , −18.9 •+16.6 −11.5 ).

Fig. 9 . 90 5DQGRPSRLQWFig. 10 .
Fig.9.Distribution of discrepancy degree D max in 500 simulated isotropic datasets.The upper two panels show the results of statistical isotropic analyses (isotropy).The lower two panels show the results of statistical isotropic analyses that preserve the spatial inhomogeneity of real data (isotropy RP).The purple and blue represent the statistical results of Ω m and H 0 .The black curve is the best fit to the Gaussian function.The solid black and vertical dashed lines are commensurate with the mean and the standard deviation, respectively.The red lines represent the discrepancy degree from the real data.For the isotropy analyses, the statistical significance of the real data are 3.68σ for Ω m anisotropy and 3.01σ for H 0 anisotropy.For the isotropy RP analyses, the statistical significance of the real data are 2.09σ and 1.35σ for Ω m and H 0 anisotropy, respectively.

Fig. 17 .
Fig. 17.Relationship between average error and discrepancy degree D max and SN number used.The upper panels show the relationship between the average error and the used SN number.The lower panels show the relationship between the discrepancy degree D max and the SN number used.Blue and red correspond to parameters Ω m and H 0, respectively.

Fig. 18 .
Fig. 18.Distribution of preferred directions (l, b) with 1σ range in other independent observations.We mark the directions given by Ω m with plus signs at different screening angles, and the directions given by H 0 are shown with triangles at different screening angles.The color represents the important result obtained in this work; black shows the result given by other independent 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).

Table 1 .
Detailed information of analysis results from the all-sky distribution of cosmological parameters.D max represents the degree of deviation from the cosmological principle; the larger the value, the higher the degree of deviation.c α indicates the statistical significance of real data from the isotropy analysis.d β indicates the statistical significance of real data from the isotropy RP analysis.For the RF (60 • ) results of the Pantheon+ sample, the reduced chi-square (χ 2 r ) corresponding to Ω m,min is 0.69, which is much smaller than 1.00.Therefore, the RF (60 • ) results are not given.
a θ max is the screening angle.b *

Table 2 .
Detailed information on the investigation of the relationship between the discrepancy degree D max and z max .
Yang et al. (2014)aropoulos (2010)statistical significance.In addition, we find that D max changes with z max and has a maximum value near z max = 0.30.The average error and D 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 byAntoniou & Perivolaropoulos (2010),Cai & Tuo (2012),Wang & Wang (2014),Yang et al. (2014),Kalus et al.