Cosmic Bulk Flow Analysis in Modified Gravity Theories: f(R) and Perturbed f(R) Models with Neutrino Coupling

The accelerating expansion of the Universe, a cornerstone discovery in modern cosmology, is attributed to the enigmatic force known as dark energy. This phenomenon was first evidenced by observations of distant Type Ia supernovae (SNe) in the late-1990s, which revealed that these cosmic explosions were farther away than predicted by a decelerating universe model, implying an accelerating expansion (A. G. Riess et al. 1998; S. Perlmutter et al. 1999). Subsequent measurements of the cosmic microwave background radiation and large-scale structure surveys have corroborated these findings, suggesting that dark energy constitutes approximately 68% of the total energy density of the Universe (Planck Collaboration 2018). The exact nature of dark energy remains one of the most profound mysteries in physics, with leading hypotheses ranging from a cosmological constant (Λ) with a constant energy density to dynamic fields such as quintessence.

The Lambda Cold Dark Matter (ΛCDM) model, while highly successful in describing the large-scale structure and evolution of the Universe, is not without its challenges and unresolved issues. One of the most significant problems is the “cosmological constant problem,” which arises from the discrepancy between theoretical predictions and observational measurements of the cosmological constant, Λ. Quantum field theory predicts a value for the vacuum energy density that is about 120 orders of magnitude larger than the observed value, leading to a severe fine-tuning problem (S. Weinberg 1989; J. Martin 2012). Another issue is the “cosmic coincidence problem,” which questions why the densities of dark energy and dark matter are of the same order of magnitude precisely at the current epoch in the history of the Universe. This seems highly unlikely and suggests a need for a deeper understanding of the mechanisms driving cosmic evolution (I. Zlatev et al. 1999).

f(R) gravity, a modification of general relativity (GR), proposes an alternative to dark energy by modifying the Einstein–Hilbert action to include a general function of the Ricci scalar, R. This theoretical framework allows for a natural explanation of the accelerated expansion of the Universe without invoking dark energy. By extending the action to , where f(R) is an arbitrary function of R, the resulting field equations introduce additional terms that can effectively mimic the effects of a cosmological constant or dynamic dark energy components. This approach not only addresses the cosmological constant problem but also provides a flexible model to fit observational data across different scales. This concept has been extensively explored and discussed in various studies (e.g., M. M. Phillips 1993; S. Capozziello 2002; S. Capozziello et al. 2003; S. Nojiri & S. D. Odintsov 2003, 2006; A. De Felice & S. Tsujikawa 2010; T. P. Sotiriou & V. Faraoni 2010; S. Capozziello & M. De Laurentis 2011).

Neutrinos, the nearly massless and elusive particles, play a significant role in cosmology, influencing the large-scale structure and the Universe’s evolution. Although often considered negligible in standard gravity models, coupling neutrinos with alternative gravity theories opens new possibilities for exploring their impact on cosmic structure formation and evolution (M. Yarahmadi & A. Salehi 2023a, 2023b). These effects can be observed through cosmological phenomena, such as galaxy formation and the broader distribution of matter. Neutrinos’ interactions with f(R) gravity can also yield unique impacts on the cosmic web, introducing deviations that shape the cosmic landscape in novel ways. Investigating the unusual velocities of galaxies, which result from gravitational interactions with surrounding matter, offers a method to differentiate between the ΛCDM and f(R) theories. These velocities, deviating from Hubble’s law, serve as records of gravitational interactions, providing empirical evidence for testing these cosmological models. The concept of “bulk flow” is a significant indicator of the unusual velocity field, reflecting the large-scale gravitational forces affecting cosmic particle movement. Determining peculiar velocities accurately requires precision and distance indicators independent of redshift, such as Type Ia SNe.

The phenomenon of bulk flow provides valuable insights into the peculiar velocity field of galaxies, averaged over a large volume, and aids in understanding oscillations in matter density on a cosmic scale. The bulk flow velocity vector, calculated with high precision within specific regions, reveals complex patterns of mass distribution that extend beyond survey boundaries. This makes bulk flow a significant tool for investigating cosmological structure formation models, as it can detect mass distribution in regions not directly surveyed (J. L. Tonry & M. Davis 1981; M. Aaronson et al. 1982; D. Lynden-Bell et al. 1988; E. J. Shaya et al. 1995; C. M. Springob et al. 2007; R. B. Tully et al. 2008). Despite the valuable information it provides, the underlying cause of bulk flow remains an enigma. One critical clue is the dipole anisotropy of the cosmic microwave background (CMB), known as the CMB dipole, which reflects the motion of our Local Group (LG). The dipole amplitude is measured at 627  ±  22 km s⁻¹, oriented toward galactic coordinates l = 276° ± 3°, b = 30° ± 3° (A. Kogut et al. 1993; E. K. Conklin 1969; P. S. Henry 1971; G. F. Smoot et al. 1977; C. H. Lineweaver 1996). Initial measurements suggested local overdensities like the Virgo cluster as the cause, but further studies pointed to larger structures such as the Shapley Supercluster (A. Yahil et al. 1980, 1986; M. Aaronson et al. 1982; M. Davis & J. Huchra 1982; M. Davis & J. Peebles 1983; J. Villumsen & M. Strauss 1987; D. Lynden-Bell et al. 1989; M. Rowan-Robinson et al. 1990; O. Lahav et al. 1990; M. A. Strauss et al. 1992; E. J. Shaya et al. 1995; M. J. Hudson et al. 2004). Kashlinsky et al. reported a significant bulk flow over d ≥ 300h⁻¹ Mpc, challenging the predictions of the ΛCDM model (M. J. Hudson et al. 2004; A. Kashlinsky et al. 2009, 2011, 2012, 2010; R. Watkins et al. 2009; G. Lavaux et al. 2010; J. Colin et al. 2011; E. Macaulay et al. 2012; U. Feindt et al. 2013). Conversely, some studies found no bulk flow evidence, aligning with Cold Dark Matter (CDM) cosmology (S. Courteau et al. 2000; A. Nusser & M. Davis 2011; S. J. Turnbull et al. 2012; Y.-Z. Ma & D. Scott 2013). Discrepancies in results may be due to different underlying cosmological models, which show degeneracy at lower redshifts and deviations at higher redshifts. This study aims to measure the bulk flow of the Universe using the f(R) gravity model by analyzing variations in the luminosity distance of Type Ia SNe. Previous research has examined bulk flow in both f(R) gravity and ΛCDM models, but this study uses unique methods and different catalogs, extending its scope beyond 50h⁻¹ Mpc with redshift tomography J. Seiler & D. Parkinson (2016).

The investigation of bulk flows in the Universe provides critical insights into the underlying dynamics of large-scale cosmic structures and the fundamental forces shaping them. Modified gravity theories, such as f(R) gravity, and their perturbed versions offer a robust framework for exploring deviations from GR and their implications for cosmic evolution. The inclusion of neutrinos in these models adds an additional layer of complexity and realism, accounting for their known influence on cosmic dynamics. Neutrinos, as fundamental particles with substantial impact on the matter-energy content of the Universe, are expected to play a significant role in shaping bulk flows. By employing f(R) gravity and its perturbed variants, we can probe how modifications to gravity and neutrino interactions influence the velocity and direction of bulk flows across different redshifts. This approach allows us to test theoretical predictions against observational data, enhancing our understanding of large-scale cosmic structures and providing a more comprehensive view of the Universe’s evolution. Moreover, aligning our models with prominent cosmic structures, such as the Sloan Great Wall (SGW) and the King Ghidorah Supercluster (KGSc), helps in constraining the parameters of modified gravity theories and neutrino effects, offering valuable insights into the interplay between dark energy, gravity, and neutrinos. This investigation not only contributes to the theoretical framework of cosmology but also aids in refining observational techniques and data interpretation in the quest to unravel the Universe’s mysteries.

f(R) gravity is another intriguing concept in the world of theoretical physics. It is an alternative theory to Einstein’s General Theory of Relativity, which is the basis of our current understanding of gravity. In f(R) gravity, the “f” stands for a function of the Ricci scalar, a mathematical expression related to the curvature of spacetime (P. G. Bergmann 1968; H. A. Buchdahl 1970; A. De Felice & S. Tsujikawa 2010; H. Liu et al. 2012). The action for an f(R) gravity model in the attendance of matter components are given by

where R is the curved scalar. The equations for motion are:

which is . For the Robertson–Walker flat metric, we have:

where R₀ represents the scalar curvature corresponding to the nonperturbation metric, ρ₀ = ρ_m + ρ_ν, f₀ = f(R₀) and prim means derivative with respect to time ratio η. Combining Equations (3) and (4) for a = 1(today), we obtain,

Finally, we come to the equation of conservation:

2.1. Scalar Perturbation

Consider a flat FRW metric scalar perturbation at the length and specific timescale:

and are scalar disorders. The disturbed components of the energy-momentum tensor in this module are as follows:

where V represents the potential value for velocity disturbances. The first-order disturbed equations, assuming the background equations are kept, are as follows:

In above relations, and = ∇_α∇^α the invariant derivative of the metric ratio is not disturbed. The first-order disturbed equations in the universe during dust dominate, are obtained as follows:

The actions for an f(R) gravity model in the attendance of matter components are given by

where R is the curved scalar. The equations for motion are:

where is the effective energy-momentum tensor that includes contributions from both the standard matter and the neutrino-matter coupling and . We introduce . The term represents the functional derivative of the interaction Lagrangian with respect to the neutrino field. In the context of cosmological studies involving modified gravity and neutrinos, the coupling term in the continuity equation plays a pivotal role in elucidating the intricate interaction between these cosmic components. In the framework of f(R) gravity, where modifications to the Einstein–Hilbert action are considered, the modified continuity equation for neutrinos can be expressed as:

Here, ρ_ν denotes the energy density of neutrinos, P_ν is their pressure, H represents the Hubble parameter, and Q signifies the strength of the coupling term. The novel addition of the coupling term, Q(u^μ∇_μf_R), encapsulates the interaction between neutrinos and the modified gravity scalar (f_R). This interaction is crucial in understanding the modified dynamics of neutrinos within the cosmological context (M. Yarahmadi et al. 2025).

In simpler terms, the equation portrays how the energy density and pressure of neutrinos are influenced by the modified gravity framework. The coupling term introduces a mechanism through which neutrinos respond to the modifications in the gravitational sector described by f(R) gravity. This nuanced perspective contributes valuable insights into addressing cosmic phenomena, such as the Hubble tension, by accounting for the joint impact of modified gravity and neutrino interactions on the cosmic evolution. Our analysis not only enriches our theoretical understanding but also provides a foundation for reconciling observational data with the predictions of modified gravity scenarios in the cosmic landscape. In this paper, we use the specific form of the external source term Q_ν = − Γρ_ν where Γ = u^μ∇_μf_R.

This form is commonly used as a parameterization to account for interactions or processes affecting neutrinos in cosmological models.

The Hu–Sawicki R + f(R) gravity model stands out as a compelling modification of GR, offering significant advantages in addressing the dark energy problem and the accelerating expansion of the Universe. One of its primary benefits is its ability to naturally reproduce the late-time cosmic acceleration without the need for a cosmological constant, thus alleviating the fine-tuning issues associated with ΛCDM models. The Hu–Sawicki model also maintains consistency with solar system tests of gravity by introducing a mechanism, the chameleon effect, that suppresses deviations from GR in high-density environments while allowing significant modifications on cosmological scales. Moreover, the model provides a framework to explore large-scale structure formation, giving rise to distinct signatures in the growth rate of cosmic perturbations and the potential for alleviating the Hubble tension. By adjusting its parameters, the model can also predict deviations in the matter power spectrum, making it a useful tool for future precision cosmology and probing deviations from standard gravity.

4.1. Motivation for Choosing the Hu–Sawicki Model

Our choice of this model is motivated by the following key features:

1. 1.  
   Consistency with Observational Constraints. The Hu–Sawicki model is designed to recover GR at high curvatures, such as in the solar system, while introducing deviations from GR on cosmological scales. This property ensures compliance with stringent local gravity tests while allowing for meaningful modifications to late-time cosmic acceleration.
2. 2.  
   Phenomenological Flexibility. The model is parameterized in a way that allows for significant freedom in adjusting the strength of modifications to GR. This flexibility enables us to explore a wide range of effects on cosmic bulk flow, CMB power spectra, and structure formation, making it particularly suitable for our study.
3. 3.  
   Well-Behaved Perturbations. The Hu–Sawicki model avoids certain theoretical issues that arise in other f(R) models, such as the presence of ghosts or instabilities in the perturbation equations. The stability of the model makes it a robust framework for analyzing linear and nonlinear perturbations, as required in our analysis of bulk flow.
4. 4.  
   Compatibility with Neutrino Coupling. The Hu–Sawicki model provides a straightforward setting to incorporate and study the effects of neutrino coupling in modified gravity. Its functional form simplifies the computation of perturbed quantities and their impact on observable cosmological phenomena.

The functional form of f(R) in the Hu–Sawicki model is given by:

where:

1. 1.  
   m² is a mass scale related to the cosmological constant, which is m² = 24H₀ and H₀ = 100h.
2. 2.  
   c₁ and c₂ are dimensionless parameters, and must be best fitted.
3. 3.  
   n is a positive integer, and is often chosen to be n = 4 for phenomenological reasons.

In the context of the Hu–Sawicki f(R) gravity model, we have determined the best-fit parameters as follows: c₁ = 0.00125, c₂ = 0.0000656, f₀ = 0.99, and n = 4.

The parameter c₁ represents a coefficient that characterizes the strength of the deviation from the Einstein–Hilbert action, contributing to the modification of the gravitational dynamics on cosmological scales. In this case, c₁ = 0.00125 indicates a relatively small but significant modification, influencing the effective dark energy density and the evolution of the Universe.

Similarly, c₂ is another coefficient that works in conjunction with c₁ to further refine the f(R) function’s form. The value c₂ = 0.0000656 demonstrates the subtle adjustment required to accurately fit observational data, ensuring that the model aligns with current constraints and cosmic observations.

The parameter f₀ = 0.988 denotes the value of the function f(R) at present

Finally, the parameter n = 4 governs the behavior of the f(R) function at high-curvature regimes. This exponent controls the functional form of f(R) and determines the strength of the modification to GR in the high-curvature regime. The chosen value n = 4 ensures that the model can accommodate a range of cosmic phenomena while remaining consistent with observational constraints.

The first derivative of f(R) with respect to R, denoted f_R, is:

The second derivative of f(R), denoted , is:

In cosmological perturbation theory, these derivatives appear in the modified field equations, influencing the growth of large-scale structures, the evolution of the Hubble parameter, and other key cosmological quantities. To modify the given equation in the context of the Hu–Sawicki f(R) gravity model, we need to account for the specific form of f_R in this model, which depends on the functional shape of f(R). In the Hu–Sawicki model, f_R evolves dynamically, and thus the coupling term Γ, which is defined as Γ = u^μ∇_μf_R, will be influenced by the derivatives of f(R).

Given this, the conservation equation for neutrinos can be modified as:

where Γ is explicitly written in terms of the Hu–Sawicki model:

Thus, the modified conservation equation becomes:

The equation is similar in structure to the standard continuity equation for neutrinos in GR, but it accounts for the modifications introduced by the f(R) gravity theory.

The complete set of equations that describes the general linear perturbations for the model have been presented in the previous section. These equations are a set of nonlinear second order differential equations with a large number of variables for which there is no analytical solution except for simplest cases and only numerical analysis can be performed. Our purpose is to convert the second order differential equation to first order by introducing some new variables. There are various reasons for doing this, one being that a first-order system is much easier to solve numerically. In addition, it allows us to investigate the behavior of the system in phase space. Phase planes are useful in visualizing the behavior of the system, particularly in oscillatory systems where the phase paths can spiral in toward zero, “spiral out” toward infinity, or reach neutrally stable situations called centers. This is a useful method to determine whether dynamics of a system are stable or not. The structure of phase space of the field equations is simplified by defining a few variables and parameters,

Using the Hu–Sawicki form of f(R), the previously defined autonomous equations are modified as follows:

1. 1.  
   For ξ₁,

   

   where ξ₁ represents the normalized derivative of the scalar field.
2. 2.  
   For ξ₂,

   

   where ξ₂ is the ratio of the curvature to the Hubble parameter.
3. 3.  
   For ξ₃,

   

   where and β is defined by:

   

   and f_R and are calculated using the Hu–Sawicki model’s form of f(R).
4. 4.  
   For ξ₄,

   

   where ξ₄ describes the evolution of perturbations in the scalar field.
5. 5.  
   For ξ₅,

   

   where ξ₅ corresponds to the matter density contribution.
6. 6.  
   For ξ₆,

   

   where ξ₆ describes the gravitational potential’s evolution.
7. 7.  
   For ξ₇,

   

   where ξ₇ is related to the gravitational potential.
8. 8.  
   For ξ₈,

   

   where ξ₈ involves the neutrino contribution and Γ is a coupling term. For perturbed f(R) gravity, we remove ξ₈ and then constrain the parameters.

Where N = lna, thus . In addition, we have used the following parameters

After some calculation from equations, for simplicity, we can obtain the above parameters in terms of the new variables as

Now, for the autonomous equations of motions, we obtain

The parameter ₁ is of great importance because it allows the expression of fundamental cosmological parameters such as the deceleration parameter q and the effective equation of state (w_eff) in terms of it. Specifically, and . The deceleration parameter is a dimensionless parameter that characterizes the rate at which the expansion of the Universe is slowing down. It is defined as the negative of the ratio of the cosmic acceleration to the cosmic expansion rate squared. Mathematically, it is expressed as: , where a is the scale factor of the Universe, represents the first derivative of the scale factor with respect to cosmic time, and a^″ represents the second derivative.

The behavior of neutrinos and elusive and enigmatic particles are significantly influenced by the modified gravitational theories encapsulated in perturbed f(R) gravity (S. Nojiri & S. D. Odintsov 2011). In this section, we study constraint on the total mass of neutrinos. The equation relates the energy density of neutrinos (ρ_ν), the scale factor (a), and the modification term (1 + f_R)is given by A. De Felice & S. Tsujikawa (2010) and T. P. Sotiriou & V. Faraoni (2010):

To estimate the mass of neutrinos using this equation, additional information or assumptions are necessary. The energy density of neutrinos is expressed in terms of their mass (m_ν) and temperature (T_ν). A common approach involves using the Fermi–Dirac distribution for relativistic neutrinos (E. W. Kolb & M. S. Turner 1990; J. Lesgourgues & S. Pastor 2006):

where g_ν denotes the number of degrees of freedom for neutrinos, with a value of 2 for each neutrino species. Substituting this expression for ρ_ν into the original equation, we obtain:

To solve for the neutrino mass (m_ν), specific values for parameters and potentially additional assumptions are needed. It is important to note that cosmological models can vary, and the approach may depend on the assumptions made in the model. The relationship between Ω_ν and the sum of neutrino masses is given by J. Lesgourgues & S. Pastor (2006):

To determine the total mass of neutrinos, ∑m_ν, we need to find the best-fitting values for the cosmological parameters Ω_ν and h. The best-fitting values of Ω_ν and h are determined by identifying the region where the likelihood is maximized. These optimal values are then used in the cosmological model to calculate the elusive total mass of neutrinos, shedding light on their contribution to the cosmic structure.

The constraints obtained for the sum of neutrino masses (∑m_ν) at the 95% confidence level are ∑m_ν < 0.142 eV, which is in broad agreement with M. Yarahmadi & A. Salehi (2023a, 2023b). Additionally, for the parameter Γ, the best-fit value is Γ = 0.6 ± 0.25, derived from the Pantheon+ catalog. This result is in good agreement with M. Yarahmadi et al. (2025).

This study aims to scrutinize the bulk flow phenomena within the framework of f(R) gravity, focusing on the analysis of density perturbations, particularly leveraging Type Ia SNe data from the Pantheon catalog. Type Ia SNe are good tools to investigate anisotropy in the universe (A. Salehi et al. 2020). The preceding section established the governing equations for anisotropic f(R) gravity through scalar perturbation theory. The coupling between neutrinos and f(R) gravity introduces modifications to the gravitational force experienced by neutrinos, affecting their trajectories and distribution. This alteration in gravitational dynamics has profound implications for the evolution of bulk flow. Specifically, the modified gravity framework influences the growth rate of structure formation and the clustering behavior of cosmic matter.

In this section, we adopt the Bonvin approach, originally introduced by C. Bonvin et al. (2006a), to quantitatively assess the magnitude and orientation of the bulk flow. This methodology involves integrating variations in luminosity distance arising from the influence of a dipole effect. The luminosity distance, denoted as d_L(z, n), is examined in relation to both redshift (z) and orientation (n), and is represented by the equation

where the directional dependence is expressed in terms of spherical harmonics, introducing observable multipoles, C_l(z), for dipole l = 1. This study aims to provide a comprehensive understanding of the bulk flow phenomena through the quantitative analysis of these luminosity distance variations.

In the above equation, the direction averaged luminosity distance and dipole are, respectively, given by

and

where e is a unit vector denoting the direction of the dipole.

To derive a formula for (more details are found in C. Bonvin et al. 2006a), we use the luminosity distance to a source emitting photons at conformal time η in an unperturbed Friedmann Universe, . The motion of the observer leads to a Doppler effect which is the dominant contribution to the dipole,

where v_Bulk is our peculiar velocity. In addition, conformal time, η, is the source redshift although is not an observable quantity. To first order

where , , and . Here, is the comoving Hubble parameter. Inserting this in Equation (21), we obtain

and therefore

It is important to emphasize that Equation (5.7) introduces two distinct methodologies for computing the value of . The first approach involves deriving the parameter by utilizing Equations (2.3)–(2.6), assuming a homogeneous and isotropic Universe. Conversely, the second approach entails computing variations in luminosity distance without considering oscillations in the parameter. In this alternative method, the determination of is achieved through perturbation Equations (4.8)–(4.22), thereby accounting for the influence of fluctuations. The objective of this study is to ascertain both the direction and magnitude of the bulk flow, and for this purpose, the second approach is employed.

The Pantheon data set comprises 1701 supernovae (SNe) distributed across a redshift range from 0.001 < z < 2.3. To align this data set with a dipole anisotropy, a series of meticulous procedures are employed.

1. 1.  
   First and foremost, the equatorial coordinates of each SN are transformed into galactic coordinates.
2. 2.  
   Subsequently, the Cartesian coordinates of unit vectors corresponding to each SN in galactic coordinates are determined. These unit vectors encapsulate the directional information essential for characterizing the spatial distribution of SNe within the galactic framework.

By executing these steps, the data set is effectively reoriented and ready for further analysis, ensuring a coherent alignment with the dipole anisotropy under consideration.

where (l_i, b_i) is the galactic coordinates of the (i)th supernova . Also, , the unit vector in direction of dipole, is given by,

where (l, b) denotes bulk flow direction in galactic coordinate. So we find that

Next, we constrain the the direction and velocity of the bulk flow across different redshift ranges,by employing χ²,

where

In this section, we engage in the analysis of distinct redshifts, each considered independently without association with others. This particular approach is referred to as redshift tomography. The utilization of redshift tomography allows for a focused examination of individual redshift interval values, enabling a more detailed and nuanced understanding of the diverse characteristics associated with each specific redshift. This approach enhances the granularity of our exploration, providing valuable insights into the unique features and phenomena associated with distinct points in the cosmic timeline. We investigate the bulk flow direction and amplitude of bulk velocity for several redshifts, such as 0.001 < z < 0.016, 0.016 < z < 0.027, and 0.035 < z < 0.055 for the local Universe and 0.1 < z < 0.2, 0.4 < z < 0.6, and 0.8 < z < 1.4 for the large-scale structures.

6.1. The Local Universe

In this section, we analyze the bulk flow velocity (V_bulk) and direction for three different models: f(R) gravity, perturbed f(R) gravity, and perturbed f(R) gravity coupled with neutrinos. The results of our redshift tomography are summarized in Tables 1, 2, and 3, for three redshift ranges: 0.001 < z < 0.016, 0.016 < z < 0.027, and 0.035 < z < 0.055. In addition, the bulk flow directions in these redshifts are plotted in Figures 1–4. Figure 4 presents the direction of bulk flow across different redshift ranges for both f(R) gravity and perturbed f(R) gravity models. The top panel shows the bulk flow direction for the redshift range 0.001 < z < 0.016. On the left-hand side, the results correspond to the perturbed f(R) gravity model, while on the right-hand side we observe the bulk flow direction for the f(R) gravity model. At these low redshifts, the bulk flow direction remains relatively consistent between both models, indicating minimal deviation in large-scale motion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The middle panels of Figure 4 focus on the bulk flow direction in the redshift range 0.016 < z < 0.027. Here, a more pronounced shift in the bulk flow direction is visible between the two models, with the perturbed f(R) model on the left-hand side showing a slightly more coherent and aligned bulk flow direction compared to the f(R) gravity model on the right-hand side. This suggests that the introduction of perturbations influences the flow direction, potentially aligning it more closely with large-scale cosmic structures.

The bottom panels of Figure 4 depict the bulk flow direction for the higher redshift range, 0.035 < z < 0.055. At these redshifts, the perturbed f(R) gravity model (left-hand side) exhibits a bulk flow direction that is significantly more aligned with the large-scale cosmic motion, especially when compared to the f(R) gravity model on the right-hand side. The results indicate that at higher redshifts, the inclusion of perturbations enhances the coherence of the bulk flow direction, bringing it closer to the expected direction associated with large cosmic structures, such as the supercluster alignment.

Figure 5 illustrates the amplitude of the bulk flow for different redshift ranges across various gravitational models. The top panels show the bulk flow amplitude for the low redshift range, 0.001 < z < 0.016, where we observe a relatively moderate bulk flow. In this range, the left-hand panels represent the results for f(R) gravity, the middle panels show the perturbed f(R) gravity, and the right-hand panels display the perturbed f(R) gravity coupled with neutrinos. It is clear that the inclusion of perturbations and neutrino couplings increases the amplitude of the bulk flow.

Zoom In Zoom Out Reset image size

The middle panels focus on the redshift range 0.016 < z < 0.027, where a notable increase in the bulk flow amplitude is evident compared to the lower redshift range. Once again, the left-hand column corresponds to the f(R) gravity model, and a gradual increase in the amplitude is seen as we move to the perturbed f(R) model and the neutrino-coupled version. The rise in amplitude reflects the stronger impact of neutrino interactions at this redshift range, contributing to more coherent large-scale motions.

Finally, the bottom panels display the amplitude of the bulk flow in the higher redshift range, 0.035 < z < 0.055. At this range, the bulk flow amplitude reaches its peak, particularly in the neutrino-coupled perturbed f(R) model (right-hand panels), indicating the critical role of neutrinos in enhancing the bulk motion at higher redshifts. Across all panels, the transition from f(R) gravity to its perturbed form and then to the model coupled with neutrinos shows a consistent amplification of the bulk flow amplitude, especially as redshift increases, underscoring the influence of perturbations and neutrinos on large-scale velocity fields in the Universe.

Our results are given in Tables 1–3. From the data presented in the tables, it is evident that the inclusion of neutrinos has a significant impact on the bulk flow velocity. For the standard f(R) gravity model, the bulk flow velocities are lower across all redshift ranges compared to the models where perturbations and neutrino couplings are considered. Specifically, in the 0.001 < z < 0.016 range, the bulk flow velocity increases from 147 km s⁻¹ in the f(R) model to 173 km s⁻¹ when neutrino couplings are included. Similar trends are observed in the higher redshift bins, where the inclusion of neutrinos results in the highest bulk flow velocities, with the 0.035 < z < 0.055 range reaching a maximum of 427 km s⁻¹.

The addition of neutrinos influences the dynamics of the cosmic large-scale structure, leading to higher velocities of the bulk flow. Neutrinos, due to their relativistic nature at early times, suppress structure formation on small scales, but they also contribute to anisotropies in the matter distribution. When these neutrino effects are incorporated into the model, they enhance the overall bulk flow velocity as seen in our results. This enhancement can be interpreted as an indication that the presence of neutrinos contributes to more coherent flows of matter on large scales.

Moreover, when neutrinos are included in the perturbed f(R) gravity model, the direction of the bulk flow becomes increasingly aligned with the direction of the Local Supercluster. In the lowest redshift bin (0.001 < z < 0.016), for example, the direction of the bulk flow for the neutrino-coupled model is given by (l, b) = (306°, −12°), which is remarkably close to the direction of the GA. This alignment suggests that the addition of neutrinos not only increases the speed of the bulk flow but also influences its direction, making it more consistent with the gravitational pull exerted by the Local Supercluster.

For the intermediate redshift range, 0.016 < z < 0.027, a similar pattern is observed in the bulk flow velocities. In the f(R) gravity model, the bulk flow velocity is 375 km s⁻¹, which increases to 421 km s⁻¹ when perturbations are considered. Upon introducing neutrino couplings, the bulk flow velocity further rises to 476 km s⁻¹. This significant increase once again underscores the role of neutrinos in enhancing the coherence of matter flows on large cosmological scales.

Additionally, the direction of the bulk flow shifts closer to the Perseus–Pisces Supercluster’s direction when neutrinos are incorporated. Positioned at a redshift of approximately z ∼ 0.016, the Perseus–Pisces Supercluster is recognized as one of the most expansive formations within the cosmic web. Encompassing an angular span of around 45°, this filamentary complex extends across the constellations Coma, Hercules, and Fornax, as extensively documented in studies such as M. Joeveer & J. Einasto (1978), S. A. Gregory et al. (1981), M. Yarahmadi et al. (2024), M. Yarahmadi & A. Salehi (2024), G. L. Chincarini et al. (1983), and D. J. Batuski & J. O. Burns (1985). The total luminosity of the PP Supercluster corresponds to a mass estimate of approximately 2 × 10¹⁶ solar masses. Despite its considerable cosmic scale, the gravitational influence of the PP Supercluster is distributed over a broader region, potentially resulting in a less pronounced impact on individual galaxies when compared to the more localized effects of the GA. The radial extent of the PP Supercluster, estimated to be approximately 32 Mpc, was determined at a mean redshift of z = 0.0174 through a study conducted by R. L. Shelton et al. (2022).

For the f(R) gravity model, the bulk flow is directed at (l, b) = (120°, −7°), while for the perturbed f(R) model with neutrino coupling, the direction shifts to (l, b) = (122°, −25°). This proximity to the supercluster direction strengthens the argument that the gravitational influence of the large-scale structures, along with neutrino effects, guides the flow of matter in these redshift ranges. The addition of neutrinos notably amplifies the bulk flow speed and brings the flow’s direction closer to regions dominated by massive structures.

In the high redshift range, 0.035 < z < 0.055, the effect of neutrinos becomes even more pronounced. The bulk flow velocity in the f(R) gravity model is 341 km s⁻¹, which increases slightly to 402 km s⁻¹ when perturbations are included. However, when neutrinos are coupled to the model, the bulk flow velocity reaches its highest value of 427 km s⁻¹. This indicates that neutrinos play a crucial role in the dynamics at higher redshifts, further boosting the bulk flow.

Furthermore, the direction of the bulk flow in the neutrino-coupled perturbed f(R) gravity model aligns almost perfectly with the Shapley Supercluster’s direction, shifting from (l, b) = (279°, 6°) in the f(R) gravity model to (l, b) = (315°, 35°) when neutrinos are included. The Shapley Supercluster, situated in the constellation Centaurus, holds a central position in the galactic coordinate system at approximately (l, b) = (311°, 32°). Located around 650 million light years away from the Milky Way, this colossal supercluster possesses a substantial collection of galaxies and galaxy clusters. With an estimated mass on the order of 10¹⁶ solar masses, the Shapley Supercluster exerts a profound gravitational influence on its surroundings, impacting the motion and dynamics of galaxies within its gravitational domain. The Shapley Supercluster holds a central role in deciphering the intricate dynamics underlying the formation and evolution of the cosmic web, which constitutes the vast large-scale structure of the Universe. Its colossal mass and gravitational influence shape the distribution of matter, guiding the convergence of galaxies along filaments and contributing to the intricate architecture of the cosmic web. By influencing the motion of galaxies and galaxy clusters within its gravitational domain, the Shapley Supercluster serves as a dynamic hub, offering insights into the cosmic flow of matter on a grand scale. Furthermore, the supercluster provides a unique vantage point for probing the distribution of dark matter, enabling observational investigations that contribute essential data for refining cosmological models and understanding the broader cosmic narrative.

This remarkable alignment at higher redshifts shows that neutrinos contribute significantly to guiding the bulk flow toward regions of higher gravitational attraction, such as the supercluster. The closer alignment at this redshift range suggests that neutrinos help synchronize the large-scale velocity field with the cosmic web structures, thereby making the bulk flow more aligned with the large overdensities in the Universe.

7.1. Beyond The Local Universe

The study of bulk flow involves examining the coordinated motion of galaxies and larger cosmic structures across vast distances. In regions with redshifts exceeding 0.1, the expansion of the Universe becomes a dominant factor, influencing the trajectories of galaxies on larger scales. Investigating bulk flow at these redshifts allows astronomers to trace the intricate dynamics of matter on cosmic proportions, providing insights into the underlying gravitational forces, cosmic web structure, and the influence of massive structures such as superclusters and filaments. This exploration enhances our understanding of the cosmic assembly, large-scale structure formation, and the interplay between dark and luminous matter. Utilizing advanced observational techniques and surveys, astronomers can unravel the complexities of cosmic flows beyond the local Universe.

The SGW, heralded as a colossal cosmic structure with a formidable span of 1.37 billion light years, not only commands attention for its sheer scale but also holds relevance in the intricate dance of large-scale cosmic dynamics. Discovered through meticulous analysis of Sloan Digital Sky Survey data, the SGW dominates the regions of Corvus, Hydra, and Centaurus, constituting a remarkable of the observable Universe’s diameter. While it stands as the sixth-largest cosmic object, its interaction with the broader cosmic landscape, especially within the redshift range 0.1 < z < 0.2, influences the observed bulk flow. The SGW’s immense mass and gravitational influence become key players in shaping the gravitational dynamics within this redshift interval, potentially steering the direction and amplitude of the cosmic flow. However, debates on whether the SGW is a chance alignment of three structures or an independent cosmic entity add complexity to our understanding of large-scale cosmic architectures, prompting ongoing research to unravel the nuanced interplay between cosmic structures such as the SGW and the observed cosmic flow within specific redshift regimes.

Large redshift surveys, exemplified by the notable 6dFGS, provide meticulous distance estimates across expansive cosmic domains, intricately enhancing our understanding of the cosmic web. Within the redshift range of 0.4 < z < 0.6, a superlative celestial structure emerges—the most massive supercluster known to date, as documented by R. Shimakawa et al. (2023). This colossal entity (aptly named the King Ghidorah Supercluster) commands attention because it resides approximately 1.3 billion light years away from Earth. The KGSc, a prodigious assembly comprising a minimum of 15 massive galaxy clusters interwoven by vast filaments, stands as a testament to the sublime intricacies of the cosmic tapestry.

The sheer magnitude of the KGSc is nothing short of awe-inspiring, boasting a staggering mass of 10¹⁶ solar masses. This extraordinary mass surpasses that of our Local Supercluster, the Laniakea Supercluster, by an order of magnitude. Extending across an immense span of about 400 Mpc (equivalent to 1.3 billion light years), the KGSc showcases the vast reach of its cosmic influence. The groundbreaking discovery of the KGSc in 2022, leveraging data from the Subaru Telescope’s Hyper Suprime-Cam (HSC) survey, has profoundly deepened our comprehension of cosmic structures.

The bulk velocity (V_bulk) results for the three different redshift tomographies, as presented in the tables, show a distinct behavior for f(R) gravity, perturbed f(R) gravity, and perturbed f(R) gravity coupled with neutrinos. In the f(R) gravity model, the bulk flow velocities remain relatively moderate across all redshift ranges, with values of 1002 ± 156 km s⁻¹ for the lowest redshift range (0.1 < z < 0.2), increasing to 2780 ± 282 km s⁻¹ for the highest redshift range (0.8 < z < 1.4). These velocities show a general trend of increasing bulk motion at higher redshifts, suggesting stronger bulk flows as we move further into the Universe.

When the perturbations to f(R) gravity are introduced, the bulk velocities slightly decrease in the lower redshift ranges, with the velocity at 0.1 < z < 0.2 dropping to 953 ± 142 km s⁻¹. However, in the higher redshift range (0.8 < z < 1.4), the velocity increases to 2840 ± 320 km s⁻¹, indicating that perturbations in f(R) gravity enhance the bulk flow motion at higher redshifts. This results are in broad agreement with Salehi et al. (2021).

The most significant change occurs when neutrinos are coupled to perturbed f(R) gravity. In this scenario, the bulk flow velocities increase across all redshift ranges. For the lowest redshift range (0.1 < z < 0.2), the velocity rises to 1125 ± 128 km s⁻¹, while in the highest redshift range (0.8 < z < 1.4) the bulk velocity peaks at 3086 ± 286 km s⁻¹. This enhancement in bulk flow is a direct consequence of the presence of neutrinos, which introduce additional anisotropic effects and lead to a more substantial bulk motion as the Universe evolves.

In the case of the perturbed f(R) gravity model coupled with neutrinos, an interesting behavior is observed. Specifically, the bulk flow velocity shows a noticeable increase compared to the noncoupled models. This can be attributed to the additional energy density provided by neutrinos, which enhances the gravitational potential, and consequently increases the amplitude of the bulk flow. The influence of neutrinos is especially pronounced at higher redshifts, where their relativistic nature plays a more significant role.

Moreover, the direction of the bulk flow in this coupled model aligns almost exactly with the direction of the superclusters, particularly at redshift ranges where prominent structures such as the SGW and KGSc are present. This alignment suggests that the inclusion of neutrinos in the cosmological model strengthens the gravitational pull toward these massive structures, resulting in a bulk flow that is directed toward the densest regions in the large-scale structure of the Universe. This precise alignment highlights the crucial role neutrinos play in shaping cosmic flows and their ability to enhance the correspondence between the observed bulk flow and the gravitational potential wells of superclusters.

Table 4 summarizes various studies of bulk flow across different redshift ranges, providing a comparison with previous findings. The results of our analysis, presented in Tables 5–7, offer new insights into the bulk flow velocity and direction at distinct redshifts within the framework of f(R) gravity, perturbed f(R) gravity, and perturbed f(R) gravity coupled with neutrinos. These results not only extend our understanding of cosmic flows but also demonstrate the impact of modifications to gravity and the inclusion of neutrinos, particularly in aligning the bulk flow direction with superclusters and increasing the flow’s amplitude.

Table 4. List of Studies of Bulk Flow

Reference	Velocity	Redshift	lo	bo	Distance
	(km s−1)		(degree)	(degree)	(h−1 Mpc)
A. Kashlinsky et al. (2008, 2009, 2010, 2011, 2012)	1000	z ≤ 0.03	287 ± 9	8 ± 6	127
R. Watkins et al. (2009)	407 ± 81	z ≤ 0.2	283 ± 14	12 ± 14	857
D. D. Kocevski & H. Ebeling (2006)	507	0.035 ≤ z ≤ 0.055	306.44	29.71	127–220
S. J. Turnbull et al. (2012)	249 ± 276	z ≤ 0.2	319 ± 18	7 ± 14	857
J. Colin et al. (2011)		0.045 < z < 0.06	287	21	168–249
R. Watkins et al. (2009)	416 ± 78	z = 0.0167	282	60	45
G. Lavaux et al. (2010)	473 ± 128	0.035 < z < 0.055	220	25	127–220
A. Nusser & M. Davis (2011)	257 ± 44	0.035 < z < 0.055	276 ± 6	10 ± 6	127–220
H. A. Feldman et al. (2010)	416 ± 78	0.015 < z < 0.06	282 ± 11	6 ± 6	40–249

Download table as:  ASCIITypeset image

Table 5. Results of Redshift Tomography for High Redshift z > 0.1 for f(R) Gravity

Redshift	Vbulk (kms−1)	l°	b°	χ2	Data Number
0.1 < z < 0.2	1002  ±  156	252  ±  17	6  ±  20	194.545213	207
0.4 < z < 0.6	2250  ±  210	355  ±  12	12  ±  16	166.439788	179
0.8 < z < 1.4	2780  ±  282	302  ±  15	− 8 ± 12	20.962546	22

Download table as:  ASCIITypeset image

Table 6. Results of Redshift Tomography for High Redshift z > 0.1 for Perturbed f(R) Gravity

Redshift	Vbulk (km s−1)	l°	b°	χ2	Data Number
0.1 < z < 0.2	953 ± 142	261 ± 21	25 ± 20	191.745213	207
0.4 < z < 0.6	2230 ± 240	359 ± 17	12 ± 10	165.229788	179
0.8 < z < 1.4	2840 ± 320	302 ± 15	− 8 ± 12	20.802546	22

Download table as:  ASCIITypeset image

Table 7. Results of Redshift Tomography for High Redshift z > 0.1 for Perturbed f(R) Gravity Coupled with Neutrinos

Redshift	Vbulk (km s−1)	l°	b°	χ2	Data Number
0.1 < z < 0.2	1125  ±  128	255  ±  22	59  ±  28	188.386547	207
0.4 < z < 0.6	2290  ±  250	332  ±  18	69  ±  18	161.548796	179
0.8 < z < 1.4	3086  ±  286	330  ±  15	− 16 ± 17	19.102154	22

Download table as:  ASCIITypeset image

Figures 6–8 demonstrate the direction of bulk flow in z > 0.1 for perturbed f(R) gravity coupled with neutrinos. In Figure 9, the bulk flow direction for perturbed f(R) gravity is shown on the left-hand side, while the results for the bulk flow direction of f(R) gravity are presented on the right-hand side. This comparison highlights how perturbations in f(R) gravity alter the flow’s direction, especially in higher redshift regions, providing deeper insights into the dynamics of cosmic structures. As we can see in Figure 8, in the redshift range 0.8 < z < 1.4, the bulk flow direction in the perturbed f(R) gravity model coupled with neutrinos points toward (l, b) = (330° ± 15°, −16° ± 17°). Remarkably, this direction coincides precisely with the observed dark energy dipole. This alignment suggests that the inclusion of neutrinos in the perturbed f(R) gravity model not only enhances the bulk flow velocity but also aligns the direction of the flow with the large-scale structure indicated by the dark energy dipole, providing a deeper connection between neutrino effects and cosmic acceleration. Figure 10 indicates that the amplitude of bulk flow in the beyond of the local Universe.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, we have investigated the properties of bulk flow across various redshift ranges in the context of f(R) gravity, perturbed f(R) gravity, and perturbed f(R) gravity coupled with neutrinos. Our analysis reveals significant insights into the dynamics of large-scale cosmic flows and their alignment with prominent cosmic structures such as the SGW and KGSc.

We observe that the inclusion of neutrinos in the perturbed f(R) gravity model leads to a notable increase in the bulk flow velocity at all redshift ranges. This enhancement is particularly prominent in higher redshift ranges, with velocities exceeding 3000 km s⁻¹ in the 0.8 < z < 1.4 range. Furthermore, the direction of the bulk flow in this model shows a striking alignment with the dark energy dipole, particularly at redshifts z > 0.4, where the bulk flow direction coincides almost exactly with the direction of the cosmic superclusters. In the redshift range 0.8 < z < 1.4, for instance, the bulk flow direction aligns precisely with (l, b) = (330° ± 15°, − 16° ± 17°), which is consistent with the direction of the dark energy dipole. This result implies a strong connection between the interaction of neutrinos and the underlying structure of the Universe, potentially shedding light on the role of neutrinos in influencing cosmic acceleration.

At lower redshifts, such as 0.1 < z < 0.2, the bulk flow is aligned with the direction of the SGW, while in the 0.4 < z < 0.6 range the flow aligns with the KGSc. This correlation between the bulk flow direction and major cosmic structures suggests that large-scale anisotropies are driven by gravitational forces associated with massive cosmic structures. The addition of neutrinos in the perturbed f(R) gravity model further sharpens this alignment, providing a more accurate reflection of observed large-scale flows in the Universe.

For redshifts below z < 0.1, we observe smaller bulk flow velocities, but the flow direction remains consistent with nearby cosmic structures. In the redshift range 0.001 < z < 0.016, for instance, the bulk flow is influenced by structures like the Local Supercluster, with velocities typically lower than at higher redshifts. Even in this low redshift range, the introduction of neutrinos in the perturbed f(R) gravity model still leads to a subtle increase in bulk flow velocity, suggesting that neutrino interactions play a role in cosmic flows, even at small scales. Additionally, the bulk flow direction for z < 0.1 in the neutrino-coupled model is aligned with the Local Supercluster and similar structures, reflecting the influence of nearby gravitational sources.

In conclusion, our results highlight the crucial role of neutrinos in shaping the velocity and direction of bulk flows in the Universe across both low and high redshift ranges. The alignment of the bulk flow with the dark energy dipole and the enhanced velocity in the neutrino-coupled model provide strong evidence of the impact of neutrinos on cosmic dynamics. These findings not only deepen our understanding of cosmic flows but also offer new avenues for probing the interactions between dark energy, neutrinos, and modified gravity models in shaping the evolution of the Universe. Future studies should further explore these connections, with the goal of uncovering the precise mechanisms through which neutrinos influence large-scale cosmic structures.

In this section, we investigate the effect of the bulk flow on CMB power spectrum. Figure 11 illustrates the CMB power spectrum (C_ℓ) as a function of the multipole moment (ℓ) for various cosmological models, including ΛCDM, f(R) gravity, perturbed f(R), and perturbed f(R) with neutrinos. The power spectrum reflects the distribution of temperature fluctuations in the CMB, and its shape is sensitive to both the underlying gravitational model and the presence of massive neutrinos. Here, we analyze the influence of these models on the CMB spectrum and their relation to the bulk flow velocity.

Zoom In Zoom Out Reset image size

A.1. Effects of Models on the CMB Power Spectrum

A.2. Bulk Flow Velocity and Its Relation to CMB Power Spectrum

A.3. Conclusion