Bulk Flow Motion Detection in the Local Universe with Pantheon+ Type Ia Supernovae

The bulk flow in the Local Universe is a collective phenomenon due to the peculiar motions of matter structures, which, instead of moving in random directions, appears to follow an approximate dipole velocity flow. We apply a directional analysis to investigate, through the Hubble-Lemaître diagram, the angular dependence of the Hubble constant H 0 of a sample of Type Ia supernovae from the Pantheon+ catalog in the Local Universe (0.015 ≤ z ≤ 0.06). We perform a directional analysis that reveals a statistically significant dipole variation of H 0, at more than 99.9% confidence level, showing that matter structures follow a dipole bulk flow motion toward (l, b) = (326.°1 ± 11.°2, 27.°8 ± 11.°2), close to the Shapley supercluster (l Shapley, b Shapley) = (311.°5, 32.°3), with velocity 132.14 ± 109.3 km s−1 at the effective distance 102.83 ± 10.2 Mpc. Interestingly, the antipodal direction of this dipole points close to the Dipole Repeller structure. Our analyses confirm that the gravitational dipole system Shapley-Dipole Repeller explains well the observed bulk flow velocity field in the Local Universe. Furthermore, we performed robustness tests that support our results. Additionally, our approach provides a measurement of the Hubble constant H 0 = 70.39 ± 1.4 km s−1 Mpc−1, at the effective distance 102.8 Mpc, z ≃ 0.025.


Introduction
Our Local Universe is manifestly inhomogeneous, plenty of matter structures, like galaxies assembled in small groups or in galaxy clusters, and cosmic voids, with various underdense matter contents and sizes, around us [Rubin, 1951, de Vaucouleurs, 1953, Gregory and Thompson, 1978, de Lapparent et al., 1986, Tully and Fisher, 1987, Coles, 1996, Courtois et al., 2013, Tully et al., 2019, Avila et al., 2019, Franco et al., 2024].Peculiar velocities of cosmic objects are the result of gravitational fields, which, in turn, are originated by the surrounding matter distribution [Peebles, 1980, Kaiser, 1987].Therefore, peculiar motions are one of the best tracers of matter density fluctuations in the universe.
Near the Local Group of galaxies, to which the Milky Way belongs, a large underdensity termed the Dipole Repeller (DR) was recently discovered [Hoffman et al., 2017].Another prominent feature in our cosmic neighborhood is the Shapley supercluster [Raychaudhury, 1989, Scaramella et al., 1989], the largest cluster of galaxies of the Local Universe.
The system Shapley-DR acts, approximately, 1 arXiv:2405.11077v2[astro-ph.CO] 31 May 2024 as a gravitational dipole system producing a bulk flow motion of matter in the Local Universe [Hoffman et al., 2017].In this scenario the DR acts as it were a repulsive body causing the evacuation of matter in its surroundings; on the other side, the Shapley supercluster acts as the dominant matter attractor, pulling matter structures from every side.One important motivation for studying this gravitational dipole system, and the features of the induced bulk flow, is due to the dynamical effects produced as large peculiar velocities in our nearby galaxies [Tully et al., 2019, Peterson et al., 2022], information needed to a better calibration of the standard candles for precise measurements of H 0 at low redshifts [Scolnic et al., 2018].
Measurements resulting from some earlier studies already suggested a large bulk flow [Rubin et al., 1976, Lynden-Bell et al., 1988].More recently, reports in the literature (see, e.g., Hong et al. [2014], Kalbouneh et al. [2023], Mc Conville and Ó Colgáin [2023], Perivolaropoulos [2023]) have studied the magnitude and direction of this bulk flow velocity with diverse approaches and investigating various cosmic tracers, where these studies confirm the dipole nature of the bulk flow, but point out some differences in the velocity magnitude and/or direction [Turnbull et al., 2012, Scrimgeour et al., 2016, Qin et al., 2018, Avila et al., 2023, Watkins et al., 2023, Whitford et al., 2023].
In this work, we perform an accurate directional analysis of the Type Ia Supernovae (SNe Ia) data from Pantheon` [Brout et al., 2022] to measure the bulk flow velocity, both in intensity and in direction.Our analyses differ from other studies in various aspects.First, our methodology follows previous studies looking for preferred directions on the sky using several cosmic tracers [Bernui et al., 2008, Marques et al., 2018, Kester et al., 2024]; second, to avoid systematics from different observables, we use only SNe Ia data from Pantheon+; this choice is a challenge due to the possibility of having few SNe in some directions of the sky; third, we measure H 0 and its uncertainty in a set of N directions, which covers the celestial sphere, performing a best fit procedure using the Pantheon+ catalog and its covariance matrix, without cosmological model assumptions; fourth, we perform consistency and robustness tests to support our analyses and results.
Our directional analysis confirms an approximate dipole behaviour of the Hubble constant, with a dipole direction close to the Shapley-DR direction.The statistical significance is obtained by comparison with a set of simulated maps, produced by randomizing the angular positions of the SNe sample in study, and repeating our directional analysis procedure.For consistency, our analyses scan the celestial sphere considering two angular resolutions, but our final results and conclusions are drawn from the best angular resolution, i.e., scanning the sky with N " 192 (N side " 4) spherical caps.
The outline of this work is the following.In Section 2 we select the SNe data, subsample of the Pantheon+ SNe Ia catalog, for our directional analysis in the Local Universe.In Section 3 we describe in detail the directional analysis approach, including the measurements and uncertainties involved.In Section 4, we present the H 0 -maps, calculate their dipole component, and determine its statistical significance.In Section 5 we present our conclusions and final remarks.We leave for the Appendix our consistency and robustness tests that support our analyses.

Data: The Pantheon+ SNe Ia Catalog
The recently released catalog of SNe Ia Pan-theon+ [Brout et al., 2022], successor of the original Pantheon sample [Scolnic et al., 2018], compiles a larger number of SNe Ia events, with more precise data, and it includes events located in host galaxies whose distances were obtained using Cepheids.The Pantheon+ catalog contains 1701 type Ia SNe, with redshifts 0.001 ď z ď 2.26, and their angular distribution is shown in the top panel of Figure 1 (in galactic coordinates, where the Milky Way is in the center of the figure. In addition to encompassing several redshift reference frames, such as the cosmic microwave background (CMB), and Heliocentric, the Pan-theon+ catalog offers a comprehensive collection of precise data.This includes the light curves, covariance matrix, distance modulus, bolometric magnitude, and their associated uncertainties [Scolnic et al., 2022].In our analyses, we shall use the following information from SN Ia data: the sky angular position, the redshift, the distance modulus, the peculiar velocities, all with their re- spective uncertainties except the last one, and the covariance matrix.Regarding the peculiar velocities' uncertainties, it is worth mentioning that they are undetermined (in the catalog appears the value 250 km s ´1 for all the SNe listed there).
Our objective in this work is to detect the bulk flow motion in the Local Universe investigating the latest SNe Ia data from Pantheon+.The choice of the redshift interval for analyses was motivated by the following: (i) the gravitational dipole in analysis is limited by the region spatially close to the Shapley supercluster, where z Shapley » 0.05; for this, we consider data with z CMB ď 0.06; (ii) galaxies very close to us are not expected to follow the Hubble flow, as overdense and underdense regions in our neighborhood affect their motion, making it dominant over the Hubble expansion, then we adopt a conservative low-redshift limit z CMB ě 0.015 [Riess et al., 2016, Peterson et al., 2022].In this case, our analyses are restricted to the Local Universe; therefore one can use the Hubble-Lemaître (HL) law to find H 0 [Visser, 2004].Notice that, for consistency, in Appendix A we consider other four redshift intervals, close to z CMB P r0.015, 0.06s, where we apply our directional analysis to study the SNe samples contained in these intervals; as observed there, we basically obtain the same results.
According to these selection criteria, the sample chosen for analyses has 501 SNe Ia with redshifts 0.015 ď z CMB ď 0.06, i.e., measured in the CMB reference frame.The angular distribution of this sample is shown in galactic coordinates in the lower panel of Figure 1, and its redshift distribution is given in Figure 2.

Methodology
In this section, we present our methodology to perform the directional analysis of the set of 501 SNe Ia selected in the previous section from the Pantheon+ catalog.Our study consists in determining the value of H 0 in different sky directions.For this we define a set of N directions, chosen according to the HEALpix scheme, that uniformly scan the celestial sphere with spherical caps of angular radius γ centered at these N directions.Our directional analysis considers two angular resolutions to scan the celestial sphere: with N " 48 (N side " 2) and N " 192 (N side " 4) spherical caps.Then, we perform the best fit analysis in the HL diagram for those SNe Ia inside the ith spherical cap, centered at direction i, obtaining the best fit value H b-f 0 Ñ H i 0 at the ith direction.We repeat this procedure for the N directions, obtaining the data set tH i 0 u, for i " 1, 2, ¨¨¨, N .The details of this analysis will be discussed in Section 3.2.1.

Finding H 0 with the Hubble-Lemaître Law
In the CMB reference frame, the measured velocity of any cosmic object is [Courteau andvan den Bergh, 1999, Baumann, 2022] where peculiar motions contribute to deviate cosmic objects from the Hubble flow (i.e., the universe expansion).Considering the radial component of this equation, one obtains where D is the radial physical distance to the cosmic object, and v r pec is the radial projection of its peculiar velocity (to simplify the notation, from now on, we write: v pec " v r pec ).The way we adopt to minimize the impact of the peculiar velocities is to choose for analysis SNe with z CMB ě 0.015, because in such a case, the peculiar velocities are a small fraction of the radial recession velocity.In fact, as observed in the plots shown in Figure 3 the peculiar velocities of our sample in analysis are, on average, lower than 5% of the recession velocity at redshifts z CMB ě 0.015.The impact of the SNe peculiar velocities will be considered later in our analyses, in Section 4.2.1.

Luminosity distance
To estimate the value of the Hubble constant, H 0 , through a best fit analysis, we use the SNe Ia luminosity distance and the redshift of each SN.In the Local Universe, luminosity distance is a good approximation for estimating proper distances to objects with small redshifts, z ă 0.2 [Visser, 2004].In effect, to perform the best fit for H 0 , we use the first-order expansion of the luminosity distance derived for the HL law, resulting from the series expansion for the physical distance Dpzq " D L pzq, where z CMB is the SN redshift in the CMB frame.
To calculate the SNe Ia luminosity distances, we use the distance modulus µ, where the fidu-cial magnitude was determined from SNe host distances [Brout et al., 2022].The distance modulus is defined as the difference between the apparent magnitude m B and the absolute magnitude M B [Sparke and Gallagher, 2007] µ " m B ´MB " 5 log 10 pD L {Mpcq `25 , (4) where D L is given in Mpc.Then, we calculate the luminosity distance as D L " 10 µ´25 5 . (5)

Uncertainties associated to D L and δH 0
To determine the 1σ uncertainty of a function X that depends on two or more variables u, v, ..., X " f pu, v, ...q, we calculate it by quadrature where the terms BX Bu , BX Bv , ..., are the partial derivatives of X with respect to u, v, ..., and σ u , σ v , ..., are the uncertainties associated with the measured quantities u, v, ..., respectively.According to this, the errors associated to measurements of the luminosity distance, D L , are where σ µ is directly obtained from the Pantheon+ catalog.To estimate the error of H 0 we take into account the uncertainties from the best fit procedure, which in turn uses the covariance matrix where σ b-f is the uncertainty in the best fit, and σ st corresponds to the statistical error calculated as σ i st " where n i is the number of SNe Ia observed in the spherical cap i, and σ dp corresponds to the uncertainty due to the dipole nature of the H 0 -map (see Section 4.1) uncertainty to be calculated using a Monte Carlo simulation in Appendix B.

The SNe H 0 -maps
To investigate in what direction our sample of 501 SNe Ia, with redshifts 0.015 ď z ď 0.06, show a larger value of H 0 , or otherwise does not show a preferred direction, we perform our directional analysis.For this, we shall scan the celestial sphere with N spherical caps, of radius γ, and consider two cases: N " 48 and N " 192 caps.
Let Ω γ i " Ωpθ i , ϕ i , γq P S 2 be a spherical cap on the celestial sphere, of γ degrees of aperture, with vertex at the ith pixel with coordinates pθ i , ϕ i q, that is, pθ i cap , ϕ i cap q " pθ i , ϕ i q, for i " 1, 2, ¨¨¨, N pix (for details of our directional analysis see Bernui et al. [2008]).Both, the number of spherical caps N and the angular coordinates of their centers tpθ i cap , ϕ i cap qu, i " 1, 2, ¨¨¨, N , are defined using the Healpy HEALPix1 pixelization algorithm [Górski et al., 2005].
The pixelization parameter N side provides the grid resolution that defines the quantity of pixels, N pix , in which the sphere is pixelized where N side " 2 k , and k is an integer number sometimes called order.This means that N side " 2 produces N " 48 cap centres, and N side " 4 produces N " 192 cap centers, which correspond to the angular resolutions we will use in our directional analysis.
As observed in the lower plot of Figure 1, the SNe are not uniformly distributed on the sky, this is because the Pantheon+ catalog is a compilation of SNe that are short-lived random events in the sky, not long-lived cosmic tracers that can be observed at any time with an astronomical survey.This implies that the number of SNe is not equal in all caps, a fact that we take into account in our error analysis, and in the choice of γ, the radius of the spherical caps.

3.2.1
From the H 0 -caps to the H 0 -map Our directional analysis consists in calculating H 0 in each one of the N caps that scan the celestial sphere, computation that is done performing the best fit of the HL diagram for the SNe observed in each cap, and assigning the H 0 value found to the centre of the cap.For example, suppose that in the i-th spherical cap, of radius γ, we observe n i SNe, then we perform the best fit of the HL diagram with these n i events finding the value H i 0 .In Appendix C, we investigate the correlation between the number of SNe and the best fit H 0 , where we conclude that the number of SNe observed in each of the spherical caps has a negligible impact in the dipolar behaviour of H 0 .The set of N values: tH i 0 u, for i " 1, ¨¨¨, N , are then assembled together into a full-sky map, hereafter the H 0 -map; attributing a color scale to these set of real numbers tH i 0 u we obtain a colored H 0 -map, as seen in the left panels of Figure 4.The next step is to calculate the dipole component of the H 0 -map, and establish its statistical significance with the help of simulated isotropic maps (see Section 4.3).
The radius of the spherical cap, γ, is a challenge: it cannot be small because the number of SNe observed in that caps would be so small that the best fit is dominated by statistical noise, and cannot be large that one loses precision2 in determining the directions where H 0 is maximum or minimum.As a lower limit for the number n i we consider the number of the SNe analyzed in 1998 [Riess et al., 1998, Perlmutter et al., 1999], that is n i " 40 SNe.After various tests, we adopt γ " 60 ˝.However, we also performed robustness tests with other cap radius to verify that the value γ " 60 ˝is not biasing our results.As shown in Appendix D, we redo our directional analysis considering caps with γ " 65 ˝and γ " 70 ˝obtaining completely similar results.

Analyses and Results
In this section, we shall present our analyses and corresponding results.We start calculating the H 0 -maps, for the two angular resolutions, with N " 48 and N " 192 spherical caps.Then, we calculate the dipole components of these maps and determine their statistical significance.

Dipole structure and the angular power spectrum
With all the data in hand, we then calculate the H 0 -maps, scanning the celestial sphere with spherical caps of radius γ " 60 ˝.Our results are displayed in Figure 4.The next step is to calculate the dipole component C 1 of the corresponding angular power spectra.Clearly, the dipole term C 1 is the larger one, C 48 1 " 1.52 for the H 0 -map 48 and C 192 1 " 1.66 for the H 0 -map 192 , indicating that the H 0 -maps show a net dipolar behavior.However, one needs to evaluate the statistical significance of this result, and this is done by comparison with H 0 -maps obtained by redoing our directional analysis procedure but randomizing the choice of the SNe.The production of these sets of isotropic H Iso 0 -maps considering 48 and 192 caps, and the statistical analyses to determine the statistical significance of the dipole components of the H 0 -map 48 and H 0 -map 192 are described in Section 4.3.
Regarding the dipole interpretation of our H 0maps, we observe in Figure 4 that this is a good approximation, but not an exact result.In fact, the color of each pixel in our H 0 -maps represents the dominance of inflows or outflows, phenomena that increases or decreases the recession speed of matter structures, respectively, relative to the Hubble flow.That is, the diversity of colors in the pixels of a H 0 -map reflects the dominance of one phenomenon or the other, or even the balance between both, effects that are quantified by the best fit H i 0 from the HL diagram of those SNe in the ith cap: reddish (bluish) pixels represent H i 0 values above (below) the mean, and greenish for values close to the mean.
Remarkably, the prevalence of each of these phenomena is determined by the large-scale distribution of clustered matter and large voids, and the Local Universe is plenty of them [Courtois et al., 2013, Hoffman et al., 2017, Tully et al., 2023, Dias et al., 2023].Large matter structures besides the Shapley supercluster, like the South Pole Wall [SPW, Pomarède et al., 2020], a huge structure at z » 0.04, produces an inflow that contributes to distort a pure dipolar pattern.On the other hand, the presence of large underdense regions besides the DR [Franco et al., 2024], like the Perseus-Pisces void (z » 0.027), is producing outflows that also leave their peculiar signature in our H 0 -maps.All these complex gravitational interactions produce inflows and outflows that affect the dynamics of the SN hosts, phenomena that are captured by our directional analysis and revealed in our H 0 -maps.The observed dipole structure of the H 0 -maps is just the consequence of a dominant gravitational dipole system, i.e., Shapley-DR [Hoffman et al., 2017].

Bulk flow velocity
Given a data set of low-redshift SNe Ia, z ! 1, one can perform a best fit analysis of the linear HL law and write where the angular parentheses mean the procedure done to obtain the best fit value H b-f 0 .Clearly, according to the left panels of Figure 4, the H 0 -maps, the value one obtains in the best fit analysis, H b-f 0 , depends on the direction in study.Analyzing these maps one determines the positive (`) and negative (´) dipole directions of the H 0 -maps (right panels of Figure 4).Consider the best fit value of the H 0 -map obtained in the positive (negative) dipole direction, that is, H 0 (H 0 ).Then, the directional analysis of the SNe data in these directions, with effective distance R, states where we interpret the excess (defect) in H 0 as being the effect of the bulk flow motion, with velocity V BF pRq at the effective distance R, on the universe expansion rate measured by the SNe.As a consequence of this, one has Therefore, In what follows we will show how to obtain the effective distance R.

Effective distance R and peculiar velocities
To calculate the effective distance R of a data sample, necessary to compute the bulk flow velocity, we use the thermal noise model of Turnbull et al. [2012], where D Li is the luminosity distance of the i-th SN, σ i is the peculiar velocity uncertainty of the i-th SN, and N SNe is the number of SNe in our sample.Unfortunately, the Pantheon`catalog does not provide uncertainties of the SNe peculiar velocities measurements3 .Therefore, for the calculation of the effective distance, R, we need to estimate the set tσ i u " tσ vpec i u.
From equation (2), the peculiar velocity in the radial direction is given by To calculate the uncertainty of v pec , σ vpec , we use equation ( 6) multiplied by a factor4 1{ ?2, obtaining5 We then compute  to the positive and negative dipole direction of the H 0 -map.Performing error propagation, the uncertainty of R, σ R , is calculated from the uncertainties in the luminosity distance, σ D L i , and in the peculiar velocity, σ i , that is, We then calculate the effective distance using equation ( 15), and its associated uncertainty using equation ( 18), obtaining R 48 " 104.67 10.31 Mpc and R 192 " 102.83 ˘10.24Mpc, for N " 48 and N " 192 spherical caps, respectively.For consistency, we have also estimated the effective distance using Monte Carlo simulations, as presented in Appendix E, obtaining values fully consistent with these ones.

Bulk flow velocity with SNe Ia
Our measurements of the directional variation of the Hubble constant were done in Section 4.1, where we have obtained a net dipolar structure in the H 0 -maps, as seen in Figure 4 4, are interpreted as being originated in the combination of motions of the SNe host galaxies due to the Hubble flow (i.e., the universe expansion) and due to peculiar motions relative to the matter structures, including the bulk flow motion.Clearly, in a universe without peculiar motions, the H 0 -maps should be mono-colour, indicative of the isotropy of H 0 , that is, the same value H 0 in any direction (except statistical fluctuations originated in the best fit procedure of the HL diagrams).According to equation ( 14), to obtain the magnitude of the bulk flow velocity, we need δH 0 .We calculate δH 0 " H 0 ´H0 , considering the values in the H 0 -maps (left panels of Figure 4) corresponding to the `{´dipole directions shown in the right panels of Figure 4 5).
In Figure 6 we can observe some interesting directions on the celestial sphere, including the dipole bulk flow motion found in our analyses.Our results, summarized in Table 1, confirm that the gravitational system Shapley-DR explains well the dipole nature of the bulk flow motion, and that matter structures in the Local Universe appear to move in the direction of the Shapley supercluster with a velocity of 132.14 ˘109.30km s ´1.6

Isotropic maps and statistical confidence analysis
We want to verify that the dipole signature found in our H 0 -map analyses is not a biased result, perhaps due to the methodology of analysis, or some artifact in the data.Our statistical confidence test considers the comparison with isotropic H 0maps, as follows.
Suppose that the number of SNe in the ith spherical cap is n i , for i " 1, ¨¨¨, N .These n i SNe were selected in the ith cap because their angular separation with respect to the center of the cap is less or equal γ degrees.In this test, instead, for the ith cap we randomly choose n i SNe from any part of the sky, where i " 1, ¨¨¨, N ; then, we redo the best fit procedure in the N HL diagrams obtaining an isotropic H 0 -map, termed H Iso 0 -map.We repeat this process 1000 times, obtaining 1000 H Iso 0 -maps.For illustrative purposes, in Figure 7 we show two of these H Iso 0maps, both for N " 48 and 192 spherical caps, i.e., H Iso-48 0 -maps and H Iso-192 0 -maps, respectively.We also calculate the angular coordinates of the dipole components from these H Iso 0 -maps, and show them collectively as histograms in Figure 8, again for both angular resolutions, that is, for N " 48 and N " 192 spherical caps.

`σ2
Iso-192 , finding that this dipole is significant at more than 99.9% confidence level.

Measuring the Hubble constant
Although our directional approach is not intended to perform a H 0 measurement, even so, one can obtain this important information from the monopole components of our H 0 -maps 48{192 . 7  This allows us to calculate the Hubble constant at the effective distance R 192 " 102.83 ˘10.24Mpc, that is, z » 0.025, obtaining H 0 " 70.39 ˘1.40 km s ´1 Mpc ´1 using the H 0 -map 192 and considering spherical caps with γ " 60 ˝.For completeness, we also calculate for other angular resolutions and spherical cap sizes used in our scrutiny; the results are summarized in Table 2.
In Appendix F we discuss the influence that overdensity and underdensity matter structures can produce on H 0 measurements.
, where C 0 is the monopole term of the H 0 -maps.

Conclusions
Cosmic matter structures grow over time due to gravitational interaction.This explains why our cosmic neighborhood, z » 0, exhibits amazingly large underdense and overdense matter struc-  C 1 " 1.66 from the original maps, respectively.As observed, in both cases, our directional analysis show that the dipole components of the H 0 -maps displayed in Figure 4 are statistically significant, at more than 99.9% confidence level.
tures [Rubin, 1951, de Vaucouleurs, 1953, Courtois et al., 2013, Hoffman et al., 2017, Tully et al., 2019, Colgáin, 2019, de Carvalho et al., 2020, Avila et al., 2021, 2022, Marques et al., 2024].In this scenario, it was recently proposed that the DR and the Shapley supercluster structures dominate the galaxy velocity field in the Local Universe, acting like a dipolar gravitational system.This scenario deserves investigation because in the Local Universe the effect of peculiar motions caused by the distribution of matter competes with the expansion of the universe making the analysis of the HL diagram difficult, affecting the measurement of the Hubble constant.
Motivated by this proposal, we investigated the subsample of 501 SNe Ia of the Pantheon+ catalog in the Local Universe, with redshifts 0.015 ď z ď 0.06, aiming to measure both the direction and magnitude of the bulk flow velocity.Our directional analysis scan the celestial sphere considering two angular resolutions: with N " 48 and N " 192 spherical caps, although our conclusions are drawn from the best angular resolution case, i.e., for N " 192 caps.We found a statistically significant dipole variation of the Hubble constant, at more than 99.9% confidence level, with H 0 values in the `{´dipole directions: H Our studies confirm that matter structures in the Local Universe are indeed following a dipolar bulk flow motion toward pl, bq " p326.˝1 ˘11.˝2, 27.˝8 ˘11.˝2q, a direction close to the Shapley supercluster (l Shapley , b Shapley q " p311.53 ˝, 32.31 ˝), with the velocity of 132.14 109.30km s ´1 at the effective distance 102.83 10.24 Mpc = 71.98 ˘7.17 Mpc h ´1, using h " 0.7 to illustrate the measurements' comparison shown in Figure 5.
Not less interesting, the antipodal direction of the dipolar bulk flow found points close to the DR structure, as illustrated in Figure 6.
As shown in Figure 5, recent studies of the CosmicFlows-4 data [CF4, Tully et al., 2023] at scales greater than 100 Mpc h ´1 have reported bulk flow velocities that are in great tension with the velocities expected in the standard cosmological model, like the measurements of Watkins et al. [2023] who found 395˘29 km s ´1 at 150 Mpc h ´1 and 427 ˘37 km s ´1 at 200 Mpc h ´1, and the measurement of Whitford et al. [2023] who found 428˘108 km s ´1 at 173 Mpc h ´1, values with very low probability of occurrence: 0.015%, 0.00015%, and 0.11%, respectively (these data are also displayed in Figure 5).As suggested by Watkins et al. [2023], these results may be an indication that the matter rest frame of the universe is not that determined from the dipole in the CMB radiation [see, e.g., Migkas et al., 2021].However, several analyses performed at moderate scales À 100 Mpc h ´1 estimate the bulk flow velocities finding that they are statistically consistent with the ΛCDM prediction, as well as the measurement of Whitford et al. [2023] at 49 Mpc h ´1 analysing the CF4 dataset, and also our measurement at the effective distance of 72 Mpc h ´1.
Concluding, our directional analysis shows that the bulk flow velocity field in the Local Uni-verse is well explained by the gravitational dipole system Shapley-DR [Hoffman et al., 2017].In addition, we are confident that our results are robust, and the methodology adopted is unbiased due to the various consistency tests done, which we present in Appendix.

A Consistency test for other redshift intervals
We have investigated the consistency of our results studying the SNe present in other redshift intervals, close to the original interval analyzed in this work.Thus, we performed our directional analysis for four SNe samples with redshift intervals as follows: z P r0.01, 0.06s, z P r0.01, 0.055s, z P r0.01, 0.065s, and z P r0.015, 0.065s; the number of SNe Ia in these intervals was 565, 552, 576, and 512, respectively.The H 0 -maps 192 obtained in each case are displayed in the left column of Figure 10.To illustrate how different these maps are with respect to the original H 0 -map 192 , analyzed for the SNe with z P r0.015, 0.06s and displayed in the left panel, second row, of Figure 4, we calculate the map difference where the index i refers to each one of the four samples mentioned above.The maps tD i u, displayed in the right column of Figure 10, show that the differences with respect to H 0 -map 192{orig are very tiny for the first two cases, and differences À 2% for the last two cases, implying that our results regarding the interval analyzed in this work, z P r0.015, 0.06s, are robust.

B Uncertainty Due to the Dipole δH 0 : Monte Carlo Simulations
The dipolar nature of δH 0 pl, bq contributes to the uncertainty in the measurement of H b-f 0 (see Section 3.1.2).To calculate this uncertainty, σ dp , due to the directional dependence observed in H b-f 0 pl, bq, we perform Monte Carlo simulations.
We first generate 1000 random samples of SNe, where their angular positions pl, bq were preserved, but their luminosity distances were estimated randomly from a Gaussian distribution, considering, as the mean and the standard deviation of such distribution, the luminosity distance of each SN D L , and its uncertainty σ D L , respectively.We then performed the same process of directional analysis considering both angular resolutions, that is, with N " 48 and N " 192 caps.Then, performing the H 0 best fit procedure in the HL diagrams, we obtain 1000 random H 0 -maps Ran-48 and 1000 random H 0 -maps Ran-192 . 8 For each one of these maps, we calculate the amplitudes tδH 0 u Ran-48/192 , the dipole components tC 1 u Ran-48/192 , and the dipole directions tpl, bqu Ran-48/192 .The corresponding histograms of the distributions obtained are shown in Figure 11 (the first and second rows correspond to N " 48, and the third and fourth rows correspond to N " 192).The standard deviation of the distributions tδH 0 u Ran-48/192 provides a measure of the uncertainty due to the dipole nature of the H 0 -map, δH 0 , where we obtain σ 48 dp " 1.26 km s ´1 Mpc ´1 and σ 192 dp " 1.24 km s ´1 Mpc ´1, respectively.Another quantity that we analyze in Figure 11 are the set of dipole components tC 1 u Ran-48/192 of the H Ran-48/192 0 -maps, where their standard deviations provide a measurement of the uncertainties in    Ran-48/192 from which one can estimate the uncertainties of quantities used in our approach.In the first row we display the amplitudes tδH 0 u Ran-48 (left panel) and the dipole component tC 1 u Ran-48  (right panel) calculated from random H 0 -maps Ran-48 , where their standard deviations provide the uncertainties σ 48 dp " 1.26 km s ´1 Mpc ´1, and σ C Ran-48 1 " 0.67, respectively (see Appendix B for details).The third row show similar calculations, but for the N " 192 case, namely, σ 192 dp " 1.24 km s ´1 Mpc ´1, and σ C Ran-192   1 " 0.35.The second and fourth rows provide the distributions of the dipole directions of the H 0 -maps Ran-48/192 , showing mean and standard deviations as expected due to the procedure used in these Monte Carlo simulations, that is, pl, bq Ran-48 " p326.˝58 ˘14.˝33, 27.˝33 ˘12.˝76q and pl, bq Ran-192 " p325.˝68 ˘6.˝76, 27.˝17 ˘6.˝1q.The vertical dashed lines in the tC 1 u Ran-48/192 plots (right panels in the first and third row) represent the dipole values C 48 1 " 1.52 and C 192 1 " 1.66 of the H 0 -maps displayed in Figure 4 (see Section 4.3 and Figure 9).
C Studying the correlation between the number of SNe and the best fit H 0 An interesting question regards the possible bias of the number of SNe observed in different directions could have an impact in the best fit analysis of the HL diagram to obtain H 0 .
For this, we perform a correlation analysis of the H 0 -map with the number-of-SNe-map (that is, the map where the color in each pixel represents the number of SNe used to construct our H 0 -map, see Figure 12), using the linear Pearson correlation coefficient P. Our analyses show (the bars mean absolute value): |P| " 0.234, for the case with N " 48 spherical caps, and |P| " 0.198, for the case with N " 192 spherical caps.According to the literature, for values of the Pearson coefficient in the interval, |P| P r0.0, 0.199s means that the correlation between the pairs of maps with the same angular resolution is very weak.
Notice that, for consistency, our analyses scan the celestial sphere considering two angular resolutions: with N " 48 and with N " 192 spherical caps.However, our conclusions are drawn from the best angular resolution case, i.e., scanning the celestial sphere with N " 192 spherical caps.Due to this, we are led to conclude that the number of SNe observed in the 192 spherical caps has a negligible impact in the -statistically significant-dipolar behavior that our directional analysis found in the H 0 -map.In this appendix, we perform robustness tests to verify our directional analyses results considering spherical caps with radius different from γ " 60 ˝, the case considered along this work, namely γ " 65 and γ " 70 ˝.Moreover, we make these tests for both angular resolutions, with 48 and 192 spherical caps.
Observing Table 3, and the H 0 -maps displayed in Figure 13, we obtain an excellent agreement with the results found along the text considering spherical caps with γ " 60 ˝.Therefore, we conclude that our results, summarized in Table 1, are robust considering directional analysis done with spherical caps of different sizes used to scan the celestial sphere.

E Consistency in the Calculation of the Effective Distance
The calculation of the effective distance is an important step in our methodology of finding the velocity of the bulk flow.Equation ( 15) is based on the thermal noise model of Turnbull et al. [2012],  that needs the information of the peculiar velocities' uncertainties, which are individually unknown in the Pantheon+ catalog, and just appears as value 250 km s ´1 for all SNe.
In Section 4.2.1, we have calculated σ vpec using Equation ( 17).However, one can adopt the Monte Carlo method to find σ vpec and the effective distance R. For this, we simulate Gaussian distributions for each SN located within the spherical caps corresponding to both the `{´dipole directions (see Table 1), considering the value v i pec (data in the Pantheon+ catalog) as the mean of the distribution and 250 km s ´1 as its standard deviation.For each distribution, one for each SN in the analysis, we randomly choose a value that we consider the uncertainty σ vpec of such SN.The set of values tσ vpec u obtained in this form is termed one realization.With this set of uncertainties, we calculate the value of the effective distance R using Equation ( 15) for each realization.
We repeat this procedure to produce a set of N realiz " 5000 realizations9 , for each angular resolution, that is, for 48 and 192 spherical caps.The set of effective distances obtained tR j u, j " 1, ¨¨¨, 5000 produces two distributions for 48 and 192 caps, respectively, that we plot as histograms in Figure 14.The median of these distributions provides the value of the effective distance R. Therefore, we obtain R 48 " 118.33 ˘33.71Mpc, and R 192 " 117.57˘33.60 Mpc, results that are in good agreement with those found in Section 4.2.1, with larger uncertainties than previous ones.F The Hubble Tension: Are Overdensities and Underdensities Biasing the H 0 Measurements?
Recent works [Giani et al., 2024, Mazurenko et al., 2024, Kenworthy et al., 2019] have investigated the possibility that overdensities (i.e., large galaxy clusters) and underdensities (i.e., large voids) affect the H 0 measurements in the Local Universe, because such structures produce large inflows and outflows, respectively.In fact, this biasing effect exists, being important quantify the impact of such phenomena in the measurements of H 0 .According to Giani et al. [2024], the large supercluster Laniakea, which hosts the Milky Way, produces a negative average expansion of " ´1.1 km s ´1 Mpc ´1, inducing a variation in the Hubble constant ∆H 0 « 0.5 km s ´1 Mpc ´1.Therefore, the inflows produced by Laniakea can increase the Hubble tension when applying this correction to the SNe Ia dataset.On the other hand, Mazurenko et al. [2024] noticed that a large supervoid could be generating outflows on scales 100 ´250 Mpc h ´1 inducing a change in the bulk flow velocities at such scales, then decreasing the tension due to the measurements reported by Watkins et al. [2023] and Whitford et al. [2023]; data that we commented at the end of Section 5. Lastly, Kenworthy et al. [2019] affirm that overdensities and underdensities can change the Hubble constant in 2.2%, not affecting the Hubble tension.
In this context, our directional analysis can didactically explain these results, and also provide a quantification of this biasing effect.The color of each pixel in our H 0 -maps (see, e.g., the left column in Figure 4) is the manifestation of the dominance of inflows or outflows, phenomena that increases or decreases the recession speed of matter structures with respect to the Hubble flow.In fact, the prevalence of these phenomena in the Local Universe is determined by the large-scale distribution of clustered matter and large voids, which can vary substantially along different directions as recently noticed by Franco et al. [2024].However, this diversity is captured by our methodology, which provides, for the ith spherical cap in analysis, the best fit H i 0 from the HL diagram of those SNe in the ith cap; since the ith cap has its vertex at the position of pixel i, it will have the color corresponding to the dominant effect in that direction: reddish (bluish) for values above (below) the mean, and greenish for values close to the mean.The relative differences of the pixel values of the H 0 -map with respect to the H 0 value shown in Table 2 are À 3% in all cases, i.e., considering different angular resolutions and cap sizes, in perfect agreement with the outcomes of Kenworthy et al. [2019].

Figure 1 :
Figure 1: Distribution of SNe Ia on the celestial sphere with a Mollweide projection in galactic coordinates.Left panel: the Pantheon`sample, showing 1701 SNe Ia.Right panel: the sample of 501 SNe Ia selected for our directional analysis, within the redshift interval 0.015 ď z CMB ď 0.06.

Figure 2 :
Figure 2: Frequency histogram of the 501 SNe Ia sample selected for our directional analysis, in different redshift intervals within the range 0.015 ď z CMB ď 0.06; the bin size is ∆z " 0.002.

Figure 4 :
Figure 4: First row: the H 0 -map 48 with resolution N side " 2 (48 spherical caps, left panel), and its dipole component map (right panel).Second row: the H 0 -map 192 with resolution N side " 4 (192 spherical caps, left panel), and its dipole component map (right panel) showing the position of the Shapley supercluster as a multiplication sign X (in fuchsia color) to evidence that the dipole direction obtained in our analyses is very close to this matter structure.Note that all these maps are in units km s ´1 Mpc ´1.

Figure 6 :
Figure 6: Mollweide projection in galactic coordinates of the structures and directions characterizing flows in the Local Universe: the Dipole Repeller, Shapley supercluster, CMB dipole, and the direction of the bulk flow we found for H 0 -map 192 .

Figure 9 :
Figure9: Histograms of the dipole components, C 1 , obtained from 1000 randomized maps, plus the vertical dashed line draw for comparison from the H 0 -map 48 (left panel) and H 0 -map 192 (right panel), with C 1 " 1.52 and C 1 " 1.66 from the original maps, respectively.As observed, in both cases, our directional analysis show that the dipole components of the H 0 -maps displayed in Figure4are statistically significant, at more than 99.9% confidence level.

Figure 10 :
Figure 10: Left column: the H 0 -maps 192 for different redshift intervals.Right column: the maps of the difference between the original H 0 -map 192 (see Figure 4) and the corresponding map to the left.

Figure 11 :
Figure11: These plots show the results of our Monte Carlo simulations performed to produce H 0mapsRan-48/192 from which one can estimate the uncertainties of quantities used in our approach.In the first row we display the amplitudes tδH 0 u Ran-48 (left panel) and the dipole component tC 1 u Ran-48 (right panel) calculated from random H 0 -mapsRan-48 , where their standard deviations provide the uncertainties σ 48 dp " 1.26 km s ´1 Mpc ´1, and σ CRan-48

Figure 12 :
Figure 12: Number of SNe maps, for 48 (left) and 192 (right) spherical caps analyses.These maps show the number of SNe analyzed in each spherical cap, as in the directional analysis procedure, by assigning a color to the corresponding number of SNe inside the 60 ˝cap.These maps are useful to study the possible correlation with the maps shown in Figure 4; our results, discussed in Appendix C, confirm that the correlation is very weak.

Figure 13 :
Figure 13: Robustness tests to confirm our directional analysis results considering different spherical cap sizes, namely: γ " 65 ˝and γ " 70 ˝, tests done for both angular resolutions, with 48 and 192 spherical caps.The plots at the top correspond to the cases with 48 caps for γ " 65 ˝(left) and γ " 70 ˝(right).The plots at the bottom correspond to the cases with 192 caps for γ " 65 ˝(left) and γ " 70 ˝(right).As observed, comparing the dipole directions of these H 0 -maps, displayed inTable 3, with those in Table 1 both results are in excellent agreement.

Figure 14 :
Figure14: Histograms of the distributions obtained using the Monte Carlo approach (see the text for details) to calculate the effective distance R, in units of Mpc, for H 0 -map 48 (left panel) and H 0 -map 192 (right panel).The median values of these distributions, that provide our measurements of R, are: R 48 " 118.33 ˘33.71Mpc, and R 192 " 117.57˘33.60 Mpc, respectively.
-map, respectively (all H 0 values are in units of km s ´1Mpc ´1).The differences among diverse directional measurements of H 0 , observed in the H 0 -maps shown in the left panels ofFigure " 69.74 ˘1.30,where all H 0 values are in units of km s ´1 Mpc ´1.The effective distances, as calculated before, are R 48 " 104.67 ˘10.31Mpc and R 192 " 102.83 ˘10.24Mpc, then finding V BF 48 " 127.81 ˘110.97 km s ´1 and V BF 192 " 132.14˘109.30km s ´1, respectively.Our results are robust, and in good agreement with the values reported in the literature (see Figure

Table 2 :
H 0 measurements for different spherical cap sizes.

Table 3 :
Robustness tests to verify our directional analysis results considering different spherical cap sizes.