Null Signal for the Cosmic Anisotropy in the Pantheon Supernovae Data

The cosmological principle assumes that the universe is homogeneous and isotropic on cosmic scales. There exist many works testing the cosmic homogeneity and/or the cosmic isotropy of the universe in the literature. In fact, some observational hints of the cosmic anisotropy have been claimed. However, we note that the paucity of the data considered in the literature might be responsible for the `found' cosmic anisotropy. So, it might disappear in a large enough sample. Very recently, the Pantheon sample consisting of 1048 type Ia supernovae (SNIa) has been released, which is the largest spectroscopically confirmed SNIa sample to date. In the present work, we test the cosmic anisotropy in the Pantheon SNIa sample by using three methods, and hence the results from different methods can be cross-checked. All the results obtained by using the hemisphere comparison (HC) method, the dipole fitting (DF) method and HEALPix suggest that no evidence for the cosmic anisotropy is found in the Pantheon SNIa sample.


I. INTRODUCTION
In modern cosmology, it is usually assumed that the universe is homogeneous and isotropic on cosmic scales [1,2]. This is the well-known cosmological principle, which plays a fundamental role. Although it is indeed a very good approximation across a vast part of the universe (see e.g. [3,4]), the cosmological principle has not yet been well proven on the scales ∼ > 1 Gpc [5]. Therefore, it is still of interest to carefully test both the homogeneity and the isotropy of the universe.
The cosmic homogeneity can be broken in the well-known Lemaître-Tolman-Bondi (LTB) void model [6] and others akin to it. In such kind of models, the cosmic acceleration could be explained without the needs of dark energy or modification to general relativity. Thus, the LTB-like models have been extensively studied in the literature. On the other hand, actually the cosmic homogeneity has been tested by using various observations. For conciseness, we refer to e.g. [7,8] and references therein for details.
In the present work, we are mainly interested in the cosmic isotropy. It can be broken in the models of Gödel universe [9], Bianchi type I ∼ IX universes [10], Finsler universe [11], and so on. In fact, some observational hints of the cosmic anisotropy have been claimed in the literature. For instance, it has been found that there is a preferred axis in the cosmic microwave background (CMB) temperature map (known as the "Axis of Evil" in the literature) [12-14, 66, 69-71]. Another kind of hints for the cosmic anisotropy comes from the distribution of type Ia supernovae (SNIa) [8,[15][16][17][18][19][20][21][22][23][24][25][26][27][28]65]. It is found that there exists a preferred direction (or more) in various SNIa datasets (e.g. Union2, Union2.1, JLA), mainly by using the hemisphere comparison (HC) method proposed in [15] and then improved by [16], as well as the dipole fitting (DF) method proposed in [17]. Also, these works have been extended to include gamma-ray bursts (GRBs) [29][30][31], and the preferred direction still exists. On the other hand, it is claimed in 1998 that the fine structure "constant" α is time-varying [32,33] (see also e.g. [34][35][36][37]). A dozen of years later, a preferred direction was also found in the ∆α/α data [38,39], which means that the fine structure "constant" α is also spatially varying. Note in [17] that the preferred direction in the ∆α/α data might be correlated with the one in the distribution of SNIa. Furthermore, similar preferred direction(s) has/have also been found in other observational data, such as rotationally supported galaxies [40,41], quasars and radio galaxies [42,43], as well as the quasar optical polarization data [44,45]. Actually, in Table I and Fig. 10 of [8], we summarized the preferred directions in various observational data mentioned above. Most of them are located in a relatively small part (about a quarter) of the north galactic hemisphere (see [8] for details). In some sense, they are in agreement with each other, and hence this fact suggests that the possible cosmic anisotropy should be taken seriously.
However, it is worth noting that the paucity of the data mentioned above might be responsible for the "found" cosmic anisotropy. For example, the numbers of SNIa in the Union2 [46], Union2.1 [47] and JLA [48] datasets are 557, 580 and 740, respectively. The number of usable GRBs (without the circularity problem) is only 59 [31,49], 79 [50], or 116 [30]. The number of usable ∆α/α data [37,38] is 293. The number of SPARC rotationally supported galaxies [51] is 175. Obviously, the numbers of the data points mentioned above are ∼ O(10 2 ), which are not enough to form statistically large samples. It is reasonable to imagine that the "found" cosmic anisotropy might be just caused by statistical fluctuations, and might disappear in a large enough sample.
Very recently, the Pantheon sample [52][53][54] consisting of 1048 SNIa has been released, which is the largest spectroscopically confirmed SNIa sample to date. In fact, it is the first sample containing ∼ O(10 3 ) SNIa. Therefore, it is very interesting to test the cosmic anisotropy in the Pantheon SNIa sample, with the hope to find something different.
In this work, we test the cosmic anisotropy in the Pantheon SNIa sample by using three methods, and hence the results from different methods can be cross-checked. The rest of this paper is organized as follows. In Sec. II, we briefly introduce the Pantheon SNIa sample. In Secs. III, IV and V, the cosmic anisotropy is tested with the HC method, the DF method, and HEALPix, respectively. In Sec. VI, some brief concluding remarks are given.
Note added: When most of our computations have been completed, we were aware of a similar work by Sun and Wang [63] (SW18 hereafter) submitted to arXiv few days ago. They also tested the cosmic anisotropy in the Pantheon SNIa sample by using both the HC method and the DF method. Our work is distinct from SW18 [63] in (at least) five aspects: (1) We found the true preferred direction with the maximum anisotropy level much higher than the one found in SW18 by using the HC method. (2) We have taken the systematic covariance matrix into account when using the Pantheon data, while SW18 only considered the statistical uncertainty and ignored the systematic covariance matrix. (3) We checked the results with the simulated isotropic data in a way different from the one of SW18. (4) We found a result different from the one of SW18, by using the DF method, mainly due to the different Markov Chain Monte Carlo (MCMC) algorithms used in our work and SW18. (5) We further tested the cosmic anisotropy with HEALPix, which has not been considered in SW18. There exist other minor differences in details. We will briefly point out the differences between our work and SW18 throughout this paper. We stress that our work is completely independent, and is complementary to SW18.

II. THE PANTHEON SNIA SAMPLE
The Pantheon sample [52][53][54] consists of 1048 SNIa. The redshift range of these 1048 SNIa is given by 0.01012 ≤ z ≤ 2.26. In Fig. 1, we show the distribution of 1048 Pantheon SNIa in the galactic coordinate system (see however [64]), while the pseudo-colors indicate the redshifts of these SNIa. It is easy to see that most of 1048 Pantheon SNIa are at low-redshift. On the other hand, more than one half of 1048 Pantheon SNIa are located in the right-bottom region of Fig. 1 (see also Sec. V for details).
Instead of using the full set of parameters in the Pantheon dataset which are extremely complicated, as a convenient alternative, one can use the Pantheon plugin [54] for CosmoMC [59] to constrain cosmological models. CosmoMC is configured to expect SNIa data formatted for a SALT2 fit [55], as in the case of the JLA sample [48]), and to perform the SALT2 fit simultaneously with other cosmological constraints. The SALT2 model is an empirical model for the band correction, that relates the actual observed B-band apparent magnitudes, m B , to the inferred bolometric apparent magnitudes, m, where M B and M are the absolute B-band magnitudes and absolute bolometric magnitudes, respectively. The theoretical distance modulus is a bolometric one (over all frequencies), which is why the band correction is required. The SALT2 model [55] is linear in the stretch parameter, x 1 , and the color parameter, c, of the SNIa light curves where both the coefficients α and β are global nuisance parameters to be determined. This leads to the corrected bolometric distance modulus The Pantheon dataset has been reduced by the BBC method [67] in which the SALT2 corrections (3) are supplemented by additional corrections, namely ∆ M being a distance correction based on the host galaxy rest mass and ∆ B a distance correction from various biases predicted from simulations. Although the BBC method is not implemented in CosmoMC, one can still use CosmoMC by accepting the corrected apparent bolometric magnitudes m = µ obs + M as determined by Scolnic et al. [52,53]. As a computational heuristic [54], one sets α = β = 0 in CosmoMC (although this is of course not true of the actual fitted values [52]), and one then supplies the corrected bolometric apparent magnitudes in place of the actual B-band magnitudes expected by the CosmoMC.
As far as tests of anisotropy of the Hubble law are concerned, this means that any signature of anisotropy that is degenerate with any of the parameters in the BBC fitting will not be probed by the analysis we adopt [75]. The theoretical distance modulus predicted by the cosmological model is defined by [1,48] µ mod = 5 log 10 d L Mpc where in which E(z) ≡ H(z)/H 0 is the dimensionless Hubble parameter, and z cmb , z hel are the CMB frame redshift and heliocentric redshift, respectively. We use the Pantheon plugin [54], with the corrected bolometric apparent magnitudes m = µ obs + M , as determined from Eq. (4) by the BBC method. This leads to where for the i-th SNIa, ∆µ i = µ obs,i − µ mod,i or ∆m i = m i − m mod,i , with M being a nuisance parameter corresponding to some combination of the absolute magnitude, M , and the Hubble constant, H 0 , while the total covariance matrix C is given by [52] Here D stat is the diagonal covariance matrix of the statistical uncertainties obtained by adding in quadrature the uncertainties associated with the BBC method according to Eq. (4) of [52], while C sys is the covariance matrix of systematic uncertainties in the BBC approach, which differs somewhat from the SALT2 approach [55] since the BBC method produces distances from the fit parameters directly. The entries along with the corrected bolometric apparent magnitudes, m i , and the CMB rest-frame and heliocentric rest-frame redshifts, z cmb,i and z hel,i , are given by the data file lcparam − full − long.txt in the Pantheon plugin [54] for CosmoMC, while C sys is given in the file sys − full − long.txt [54]. Since H 0 is absorbed into M in the analytic marginalization, the Pantheon sample is independent of the normalization of the Hubble constant. According to the setting file full − long.dataset in the Pantheon plugin [54] for CosmoMC, the nuisance parameter M is marginalized by using Eq. (C1) in Appendix C of [55] (note that in the JLA dataset M is marginalized by using Eq. (C2) of [55] instead [56,57]). The best-fit model parameters (and their uncertainties) can be obtained by minimizing χ 2 Pan .
Note added: In SW18 [63], they only considered the statistical uncertainties, and actually ignored the systematic covariance matrix C sys . In the present work, we instead use the full covariance matrix with the systematic uncertainties determined by the BBC method [67], rather than simply the statistical uncertainties as in SW18.

III. TESTING THE COSMIC ANISOTROPY WITH THE HC METHOD
In the present work, we test the cosmic anisotropy in the Pantheon SNIa sample by using three methods, and hence the results from different methods can be cross-checked. For simplicity, we consider the spatially flat ΛCDM model throughout this work. As is well known, in this model, the dimensionless Hubble parameter is given by where Ω m0 is the fractional density of the pressureless matter. At first, we consider the HC method proposed in [15] and then improved by [16]. As is well known, the deceleration parameter q 0 = −1 + 3Ω m0 /2 in the flat ΛCDM model. So, it is convenient to use Ω m0 instead of q 0 to characterize the cosmic acceleration [16]. Following [16] (and e.g. [8,18,20,21]), the main steps to implement the HC method are (i) Generate a random directionr rnd indicated by (l, b) with a uniform probability distribution, where l ∈ [ 0 • , 360 • ) and b ∈ [−90 • , +90 • ] are the longitude and the latitude in the galactic coordinate system, respectively. (ii) Divide the Pantheon dataset into two subsets according to the sign of the inner productr rnd ·r dat , wherer dat is a unit vector describing the direction of each SNIa in the Pantheon dataset. So, one subset corresponds to the hemisphere in the direction of the random vector (defined as "up"), while the other subset corresponds to the opposite hemisphere (defined as "down"). Noting that the position of each SNIa in the Pantheon sample [58] is given by right ascension (ra) and declination (dec) in degree (equatorial coordinate system, J2000), one should convertr rnd andr dat to Cartesian coordinates in this step. (iii) Find the best-fit values on Ω m0 in each hemisphere (Ω m0,u and Ω m0,d ), and then obtain the so-called anisotropy level (AL) quantified through the normalized difference [16], namely (iv) Repeat for N random directionsr rnd and find the maximum AL, as well as the corresponding direction of maximum anisotropy. (v) Obtain the 1σ uncertainty σ AL associated with the maximum AL [16], Note in [16] that σ AL is the error due to the uncertainties of the SNIa distance moduli propagated to the best-fit Ω m0 on each hemisphere and thus to AL. One can identify all the test axes corresponding to AL = AL max ± σ AL . These axes cover an angular region corresponding to the 1σ range of the maximum anisotropy direction. We refer to [16] for more details of the HC method. Let us implement the HC method to the Pantheon SNIa sample [52][53][54]58]. First, we repeat 15000 random directions (l, b) across the whole sky, and find that the directions with the largest ALs concentrate around two directions (137 • , −7 • ) and (110 • , −20 • ). Then, we densely repeat 30000 random directions around these two preliminary directions. Finally, we find that the 1σ angular region with the maximum AL is in the direction and the corresponding maximum AL (with 1σ uncertainty) is In addition, we also find a sub-maximum AL in the direction (with 1σ uncertainty) and the corresponding sub-maximum AL (with 1σ uncertainty) is In fact, it is not so rare to find two preferred directions (see e.g. [8,40]). We present the pseudo-color map of AL(l, b) in Fig. 2. It is clear to see these two preferred directions within the red regions. However, it is worth noting that these two directions are not in agreement with most of the preferred directions found in various observational datasets (see Table I and Fig. 10 of [8]). This unusual fact makes these results impeachable. So, it is of interest to see whether the maximum ALs in Eqs. (15) and (17) for the real Pantheon data are consistent with statistical isotropy. To this end, we can compare the real data with the simulated isotropic data following [16]. The idea to construct the simulated isotropic data is simple. First, we fit the flat ΛCDM model to the real Pantheon data, and get the best-fit model parameter Ω m0 = 0.298111, while the corresponding M = 23.808891 (see Appendix C of [55] for technical details). Then, except that the value of m i should be simulated, we keep all the other numerical observed data unchanged for each SNIa in the Pantheon sample. We obtain m mod,i by using Eqs. (8) and (6), (11) with the best-fit Ω m0 and M mentioned above for each SNIa. Finally, the simulated m i can be obtained by taking a random number from a Gaussian distribution with the mean at the corresponding m mod,i , while the standard deviation of this Gaussian distribution is equal to the observed σ mi . So far, a simulated isotropic dataset consisting of 1048 SNIa is ready. We can then compare the absolute maximum AL of the simulated isotropic dataset with the one of the real Pantheon dataset, which are both found by repeating 1000 random directions in the corresponding dataset. We generate 100 simulated isotropic datasets and make such kind of comparisons for 100 times. In fact, 100 × 1000 computations consume a lot of computing power and time, while they have acceptable statistics. Finally, we find that the absolute maximum AL of the real Pantheon dataset is larger (smaller) than the one of the simulated isotropic dataset for 55 (45) times, respectively. This is not statistically significant in fact. Thus, the maximum AL of the Pantheon dataset is consistent with statistical isotropy. No evidence for the cosmic anisotropy is found.
Note added: In SW18 [63], they used only 500 random directions to search the maximum AL in the Pantheon sample, and hence the "maximum" AL they found is only 0.105, which is much smaller than 0.3088 found in the present work by using 45000 random directions at the price of great computing power and time. On the other hand, we checked the results with the simulated isotropic data in a way slightly different from the one of SW18 [63].

IV. TESTING THE COSMIC ANISOTROPY WITH THE DF METHOD
Let us use another method to test the cosmic anisotropy in the Pantheon data, and cross-check the result with the one obtained in the previous section. The DF method was proposed in [17], and it has been extensively considered in the literature (e.g. [8, 17, 20-22, 24, 30, 40, 41]). As mentioned in Sec. II, the observational quantity under consideration is the corrected bolometric magnitude m in the Pantheon plugin [54]. If it is anisotropic, we can consider a dipole correction, and replace m mod with wherem mod is the value predicted by the isotropic theoretical model given by Eq. (8), and A D is the dipole magnitude,n is the dipole direction,p is the unit 3-vector pointing towards the data point. In terms of the galactic coordinates (l, b), the dipole direction is given bŷ n = cos(b) cos(l)î + cos(b) sin(l)ĵ + sin(b)k , whereî,ĵ,k are the unit vectors along the axes of Cartesian coordinates system. The position of the i-th data point with the galactic coordinates (l i , b i ) is given bŷ One can find the best-fit dipole direction (l, b) and the dipole magnitude A D as well as the other model parameters by minimizing the corresponding χ 2 . In doing this, the Markov Chain Monte Carlo (MCMC) code CosmoMC [59] is used, and the nuisance parameters M can be marginalized [48]. Note that the Pantheon plugin for CosmoMC is available at [54]. Here, we let A D be a completely free parameter. In Fig. 3, we show the marginalized probability distributions of Ω m0 , the dipole magnitude A D and the dipole direction (l, b). The constraints with 1σ uncertainties are given by It is clear to see that A D = 0 is fully consistent with the Pantheon data, and the 1σ region of l and b is the whole sky. No evidence for the cosmic anisotropy is found. One can generalize Eq. (18) by including the monopole term B in addition, namely Again, it is clear to see that A D = 0 and B = 0 are fully consistent with the Pantheon data, and the 1σ region of l and b is the whole sky. No evidence for the cosmic anisotropy is found. So, by using the DF method, we confirm the result obtained by using the HC method in the previous section, namely there is no cosmic anisotropy in the Pantheon data.
Note added: In SW18 [63], they used a fairly different MCMC algorithm by their own, and obtained fairly different results. In particular, they found a relatively small 1σ region of the dipole direction, as was shown in Fig. 4 of SW18 [63]. Instead, we used the MCMC code CosmoMC [54,59], and found that the 1σ region of the dipole direction is the whole sky.

V. TESTING THE COSMIC ANISOTROPY WITH HEALPIX
Last, we consider the third method to test the cosmic anisotropy in the Pantheon data. In the HC method, the directions are compared in the way of "up hemisphere" versus "down hemisphere". So, the fine structure might be smoothed in fact. Actually, in Sec. 2 of [16] (one of the original papers proposed The whole sky is divided into 12 equal-area regions by using HEALPix, which are labeled from 0 to 11. Note that here we use the same galactic coordinate system as in Fig. 1  the HC method), the authors mentioned this issue and referred to HEALPix [60] as a solution. In e.g. [19], the cosmic anisotropy in the Union2 SNIa dataset was tested with HEALPix. Following [19], here we use it to the Pantheon SNIa dataset. HEALPix [60] is a genuinely curvilinear partition of the sphere into exactly equal area quadrilaterals of varying shape. The base-resolution comprises 12 pixels in three rings around the poles and equator. The resolution of the grid is expressed by the parameter N side which defines the number of divisions along the side of a base-resolution pixel that is needed to reach a desired high-resolution partition [60]. Following [19], we adopt N side = 1, namely the lowest resolution (12 pixels), due to the fact that the number of SNIa in the dataset under consideration is not large enough. In this case, the whole sky is divided into 12 equal-area regions by using HEALPix, which are labeled from 0 to 11, as is shown in Fig. 5. Note that here we use the same galactic coordinate system as in Fig. 1, rather than the default one of HEALPix. This can be done by using a suitable coordinate transformation.
The distribution of 1048 Pantheon SNIa in the galactic coordinate system is shown in Fig. 1. Using the routine ang2pix * provided in HEALPix package, we can find the corresponding SNIa in each region. The numbers of SNIa in all the 12 regions are given in Table I. Obviously, the numbers of SNIa in the regions 0, 4, 5, 7, 8 are N SN = 16, 0, 2, 7, 1, respectively. There are too few SNIa in these five regions. Therefore, we certainly exclude them from the following discussions. Note that the region 6 is subtle. There are 34 SNIa in the region 6. Not too many, not too few. Thus, we will consider both cases with and without the region 6. We can constrain the model parameter Ω m0 and the derived deceleration parameter q 0 = −1 + 3Ω m0 /2 by fitting the flat ΛCDM model to the corresponding SNIa in each region. In Table I, we also present the best-fit Ω m0 , q 0 with 1σ uncertainties for each region. Following [19], we use q 0 as the diagnostic for the anisotropy of the cosmic acceleration.
At first, we take the region 6 into consideration. In this case, we mask the five regions 0, 4, 5, 7, 8 mentioned above (and flag them as bad pixels). In the left panel of Fig. 6, we show the best-fit q 0 in each region, and flag the masked regions by q 0 ∼ 0. It is convenient to expand the 2D anisotropic q 0 map in spherical harmonics, namely to consider the multipole expansion [60,61], where Y m l (θ, φ) are the standard spherical harmonics, and a m l are constant coefficients which depend on the function. The term a 0 0 represents the monopole, the terms a −1 1 , a 0 1 , a 1 1 represent the dipole, and so on. Equivalently, this series is also frequently written as [61,62] q 0 (θ, φ) = a 0 + a i n i + a ij n i n j + a ijk n i n j n k + . . . , where n i represent the components of a unit vector in the direction given by the angles θ and φ, and indices are implicitly summed. The term a 0 is the monopole, a i is a set of three numbers representing the dipole, and so on. In this work, we only consider the lowest multipoles, namely the monopole with l = 0 and the dipole with l = 1, due to the low resolution of the anisotropic map. Using the routine remove − dipole * provided in HEALPix package, one can fit and remove the monopole and the dipole from a HEALPix map. Note that the masked (and bad) pixels will not be used for fit. First, we try to remove only the monopole from the left panel of Fig. 6. The best-fit monopole is −0.5756, which is equivalent to the average deceleration parameter in the whole sky. If we subtract this best-fit monopole from the q 0 map, as is shown in the middle panel of Fig. 6, the residual becomes modest (and much smaller than the monopole). Then, we try to remove both the monopole and the dipole from the left panel of Fig. 6. The corresponding best-fit monopole is −0.5316, and the three components of the dipole in the HEALPix default Cartesian coordinate system are 0.1189, 6.0550 × 10 −2 , −2.9350 × 10 −2 . This dipole is much smaller than the monopole in fact. If we subtract both the best-fit monopole and dipole from the q 0 map, as is shown in the right panel of Fig. 6, the residual becomes smaller but keeps the same magnitude as in the middle panel. The change between the right and middle panels of Fig. 6 is not so significant. The main component in the q 0 map is contributed by the monopole, and the anisotropic component contributed by the dipole is relatively small. As mentioned above, the region 6 is subtle, and it contains only 34 SNIa. Due to the paucity of SNIa in this region, the uncertainty of q 0 is fairly large, as is shown in Table I. Since we have taken the region 6 into consideration above, let us turn to the case excluding the region 6. In this case, we mask six regions 0, 4, 5, 6, 7, 8 (and flag them as bad pixels). In the left panel of Fig. 7, we show the best-fit q 0 in each region, and flag the masked regions by q 0 ∼ 0. First, we try to remove only the monopole from the left panel of Fig. 7. The best-fit monopole is −0.5440. If we subtract this best-fit monopole from the q 0 map, as is shown in the middle panel of Fig. 7, the residual becomes very small. This suggests that the q 0 map is highly isotropic in fact. Then, we try to remove both the monopole and the dipole from the left panel of Fig. 7. The corresponding best-fit monopole is −0.5301, and the three components of the dipole in the HEALPix default Cartesian coordinate system are 4.2964 × 10 −2 , 3.6163 × 10 −2 , −2.9350 × 10 −2 . The subtracted map is shown in the right panel of Fig. 7, which is very close to the middle panel of Fig. 7. In fact, the anisotropic component contributed by the dipole is very small.
Briefly, we find that there is no evidence for the cosmic anisotropy in the q 0 map. This result obtained by using HEALPix is fully in agreement with the ones obtained by using both the HC method and the DF method in the previous sections. Note added: In SW18 [63], they have not tested the cosmic anisotropy by using HEALPix. On the other hand, after the present work has been submitted to arXiv and journal, another work [68] appeared, in which the possible cosmic anisotropy in the Pantheon sample was also studied by using HEALPix. In [68], their results confirm that there is null evidence against the cosmological principle.

VI. CONCLUDING REMARKS
The cosmological principle assumes that the universe is homogeneous and isotropic on cosmic scales. There exist many works testing the cosmic homogeneity and/or the cosmic isotropy of the universe in the literature. In fact, some observational hints of the cosmic anisotropy have been claimed. However, we note that the paucity of the data considered in the literature might be responsible for the "found" cosmic anisotropy. So, it might disappear in a large enough sample. Very recently, the Pantheon sample consisting of 1048 SNIa has been released, which is the largest spectroscopically confirmed SNIa sample to date. In the present work, we test the cosmic anisotropy in the Pantheon SNIa sample by using three methods, and hence the results from different methods can be cross-checked. All the results obtained by using the HC method, the DF method and HEALPix suggest that no evidence for the cosmic anisotropy is found in the Pantheon SNIa sample.
Some remarks are in order. In this work, we only consider the spatially flat ΛCDM model. In fact, one can generalize our discussions to other cosmological models, such as wCDM, CPL models, or modelindependent parameterizations like cosmography. It is reasonable to expect that our results do not change significantly in these generalized cases.
Both the DF method and HEALPix involve the dipole and the monopole. However, in the case of the DF method, the quantity under consideration is m (which is equivalent to the distance modulus or the luminosity distance). In the case of HEALPix, the quantity under consideration is the deceleration parameter q 0 instead.
In almost all of the relevant works, only the dipole and the monopole have been considered. In fact, one can further take the quadrupole into account. Although it is not easy to deal with the quadrupole, this issue deserves future consideration.
In the case of HEALPix, we divide the whole sky into only 12 equal-area regions. This corresponds to N side = 1, namely the lowest resolution. In fact, the total number of pixels is equal to N pix = 12N 2 side . Although the number of SNIa in the dataset under consideration might be not large enough, it is still interesting to consider a higher resolution (e.g. N side = 2) in the relevant works.
One should be aware of the caveat that our conclusions might change if the anisotropic Hubble law in Eq. (18) was considered at the same time as the standardization of SNIa by the BBC method [67]. As is mentioned in Sec. II, the Pantheon sample actually use the complicated BBC method [67] (rather than the SALT2 method) in the SNIa band correction, as discussed in Sec. 3.5 of [52]. If one were to apply the SALT2 method in the manner assumed in the JLA sample with non-zero α and β, then one could apply the anisotropic Hubble law at the same time as fitting for these parameters. However, the BBC method is simply so complicated that it is not an easy task.
Nonetheless, from Fig. 1, it is easy to see that the low-redshift SNIa at z < 0.2 are relatively isotropically distributed, but the SNIa at 0.2 < z < 0.4 are concentrated in the southern latitudes with approximately  40 • < l < 190 • . Since the BBC method is designed to deal with some redshift dependent biases among other effects, empirically this could lead to degeneracies with the effects that we are searching for if the redshift ranges of the sample show an anisotropy. So, we consider a redshift tomography of the Pantheon data to explore the possible redshift anisotropy. We use the DF method in Sec. IV to find the possible dipole direction in the redshift ranges 0 − 0.2, 0 − 0.4, 0 − 0.6 and 0 − 2.26. The corresponding results are given in Table II. In all the four redshift ranges, no anisotropy has been found, since A D = 0 is fully consistent with the SNIa in each redshift range, and the 1σ region of l and b is the whole sky. This suggests that there is no (considerable) redshift anisotropy in the Pantheon sample. As is mentioned in Sec. I, there are also anomalous anisotropies in the CMB map, which still persist in the Planck data [69][70][71]. Since the CMB is extremely well sampled at all angles by the Planck satellite, these anomalous anisotropies cannot arise from lack of sky coverage as in the case of the SNIa samples. There are two possibilities that would render the CMB results compatible with the null results obtained in the present work, namely (1) the CMB results are due to anisotropies at very high redshifts, such as an intrinsic dipole on the last scattering surface [72]; (2) the CMB results are due to anisotropies from structures at very low redshifts z < 0.03 which do not give rise to a purely kinematic CMB dipole [73]. In fact, evidence for a non-kinematic dipole has been seen at the 99.5% confidence level in number counts of radio galaxies [74], which is consistent with the second possibility and the result of [43].