*Hu–Sawicki Model*

We focus on the Hu–Sawicki model with *n* = 1, considering a late-time gravity modification, described by the following function [382]:

$$f(\mathbf{R}) = \mathbf{R} + F(\mathbf{R}) = \mathbf{R} - m^2 \frac{c\_1 \left(\mathbf{R}/m^2\right)^n}{c\_2 \left(\mathbf{R}/m^2\right)^n + 1},\tag{40}$$

corresponding to the potential *V*(*φ*) in Equation (18). The parameters *c*<sup>1</sup> and *c*<sup>2</sup> are fixed by the following conditions [382]

$$\frac{c\_1}{c\_2} \approx \epsilon \frac{\Omega\_{0\Lambda}}{\Omega\_{0m}}\tag{41}$$

$$F\_{R0} \approx -\frac{c\_1}{c\_2^2} \left(\frac{12}{\Omega\_{0m}} - 9\right)^{-2},\tag{42}$$

where *FR*<sup>0</sup> is the value of the field *F<sup>R</sup>* ≡ *dF*/*dR* at the present time, and *F*(*R*) is the deviation from the Einstein–Hilbert Lagrangian density. Cosmological constraints provide <sup>|</sup>*FR*0| ≤ <sup>10</sup>−<sup>7</sup> from gravitational lensing and <sup>|</sup>*FR*0| ≤ <sup>10</sup>−<sup>3</sup> from Solar system [387,388]. We explore several choices of *FR*0.

To simplify the numerical integration of the modified luminosity distance (39), we approximate the numerical solution *yH*, obtained from the system (35), by a polynomial of order 8. This function is an accurate representation of *y<sup>H</sup>* when we restrict the solution to the range of PS (see Figure 4).

As a consequence, we obtain constraints on *c*<sup>1</sup> and *c*2, according to Equations (41) and (42). Then, we perform the same binned analysis of Section 4 using the Hu–Sawicki model and the modified luminosity distance (39).

**Figure 4.** The numerical solution for Equation (35) (blue dashed curve) plotted together with its polynomial fitting (green continuous curve) in the case of *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>7</sup> . The assumption of a function of redshift in the form of a order-8 polynomial allows an accurate fit for the numerical values. The same fitting procedure has been used in the *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>4</sup> case.

We here run the analysis both for the case of Ω0*<sup>m</sup>* fixed to a fiducial value of 0.298 and for several values of *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>7</sup> , <sup>−</sup>10−<sup>6</sup> , <sup>−</sup>10−<sup>5</sup> , <sup>−</sup>10−<sup>4</sup> (see Table 2 and Figure 5) or we let <sup>Ω</sup>0*<sup>m</sup>* vary with the two values of *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>7</sup> , <sup>−</sup>10−<sup>4</sup> (see Table 3 and Figure 6) for the SNe alone and with SNe +BAOs. Note also that the *η* parameters are all consistent for the several values of *FR*<sup>0</sup> in 1 *σ*, as you can see in Table 2, for both SNe Ia and SNe Ia + BAOs. Moreover, the values of *η* are consistent in 1 *σ* with the ones obtained from the analysis of the <sup>Λ</sup>CDM model (see also Table 1). We consider the cases *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>7</sup> , *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>4</sup> and, to study how these results may vary according to the different values of Ω0*<sup>m</sup>* chosen, we tested the model with four values of Ω0*<sup>m</sup>* = (0.301, 0.303, 0.305) taken from the 1 *σ* from a Gaussian distribution centred around the most probable value of 0.298, see [368].

We show in Figure 6 the comparison between the different applications of the Hu– Sawicki model: in the left panels (upper and lower), we consider SNe Ia only, while in the right panels (upper and lower) we combine SNe Ia+BAOs. We here remind that the assumed values for <sup>|</sup>*FR*0<sup>|</sup> of 10−<sup>4</sup> and 10−<sup>7</sup> are well constrained by the *<sup>f</sup>*(*T*) theories [167,382].

**Table 2.** Fitting parameters of *H*0(*z*) for three bins within the Hu–Sawicki model, with SNe only and SNe + BAOs with a fixed value of Ω0*<sup>m</sup>* = 0.298 and with several values of *FR*<sup>0</sup> : <sup>−</sup>10−<sup>4</sup> , <sup>−</sup>10−<sup>5</sup> , <sup>−</sup>10−<sup>6</sup> , <sup>−</sup>10−<sup>7</sup> . The columns contains: (1) the number of bins; (2) H0, (3) is *η*, according to Equation (9); (4) how many *σ*s *η* is compatible with zero (namely, the ratio *η*/*ση*); (5) *FR*<sup>0</sup> values; (6) the sample used.


Thus, the existence of this trend is, once again, confirmed, and it remains unexplained also in the modified gravity scenario. Indeed, a suitable modified gravity model which would be able to predict the observed trend of *H*0, would allow observing a flat profile of *H*0(*z*) after a binned analysis. Further analysis must be carried out with other Dark Energy models or other modified gravity theories to investigate this issue in the future, for instance focusing on the proposed model in Section 6.2.2.

**Table 3.** Fitting parameters of *H*0(*z*) for three bins within the Hu–Sawicki model, with SNe and SNe + BAOs by fixing several values of <sup>Ω</sup>*<sup>M</sup>* <sup>=</sup> 0.298, 0.303, 0.301, 0.305 and values of *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>4</sup> and *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>7</sup> . The columns are as follows: (1) the Ω0*<sup>m</sup>* value; (2) H0, (3) *η*, according to Equation (9); (4) how many *σ*s the evolutionary parameter *η* is compatible with zero (namely, *η*/*ση*); (5) *FR*0; (6) the sample used.


**Figure 5.** *Cont*.

**Figure 5.** The first four panels deal with *H*<sup>0</sup> vs. *z* for SNe, the four bottom panels include BAO measurements for the H-S model. The upper 4 panels show from the left to the right *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>7</sup> , <sup>−</sup>10−<sup>6</sup> , <sup>−</sup>10−<sup>5</sup> , <sup>−</sup>10−<sup>4</sup> , respectively. The standard ΛCDM cosmology is shown in red and the Hu–Sawicki model in blue. Analogously, the bottom panels have the same notation about the values of *FR*0.

**Figure 6.** The Hubble constant versus redshift plots for the three bins of SNe Ia only, considering the Hu–Sawicki model. **Upper left panel.** The condition of *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>7</sup> is applied to the case of SNe only, with the different values of Ω0*<sup>m</sup>* = 0.301, 0.303, 0.305. **Upper right panel.** The same of the upper left, but with the contribution of BAOs. **Lower left panel.** The SNe only case with the *<sup>F</sup>R*<sup>0</sup> <sup>=</sup> <sup>−</sup>10−<sup>4</sup> condition, considering the different values of Ω0*<sup>m</sup>* = 0.301, 0.303, 0.305. **Lower right panel.** The same as the lower left, but with the contribution of BAOs. The orange color refers to Ω0*<sup>m</sup>* = 0.301, the red to Ω0*<sup>m</sup>* = 0.303, the magenta to Ω0*<sup>m</sup>* = 0.305, and the blue to Ω0*<sup>m</sup>* = 0.298.

#### **8. Requirements for a Suitable f(R) Model**

Since the Hu–Sawicki model seems to be inadequate to account for the observed phenomenon of the decaying *H*0(*z*), in what follows, we provide some general properties that an *f*(*R*) model in the JF must possess to induce the necessary scenario of a slowly varying Einstein constant. Now, we consider again the dynamical impact of the scalar field *φ*, related to the *f*(*R*) function. Let us observe that the following relation holds in the following way:

$$\frac{d\phi}{dz} = -\frac{1}{1+z} \frac{\phi}{H} \,. \tag{43}$$

In order to get the desired behavior *φ* ' (1 + *z*) 2*η* , we must deal with a dynamical regime where the following request is satisfied:

$$\frac{\phi}{H} = -2\eta\phi\,. \tag{44}$$

We consider a slow-rolling evolution of the scalar field *φ* in the late universe, near enough to *φ* ' 1. Then, we consider in Equation (19) *ρ* ∼ 0, because we are in the Dark Energy

dominated phase, and we consider *H*<sup>0</sup> *φ*˙ small with respect to the potential term *V*(*φ* ' 1). We neglect, also, the term *φ*¨. Under these conditions, Equations (19) and (21) become

$$H^2 = \frac{V}{6\phi} \tag{45}$$

and

$$\frac{\phi}{H} = \frac{1}{9H^2} \left( 2V - \phi \frac{dV}{d\phi} \right),\tag{46}$$

respectively.

Referring to Equation (45) at *<sup>z</sup>* <sup>∼</sup> 0, we make the identification *<sup>H</sup>*<sup>2</sup> <sup>0</sup> ≡ *V*(*φ* ' 1)/6. Hence, in order to reproduce Equation (44), we must require that for *φ* → 1, the following relation holds:

$$
\eta = \frac{1}{\Im V} \left( \phi \frac{dV}{d\phi} - \mathcal{D}V \right). \tag{47}
$$

The analysis above states the general features that a *f*(*R*) model in the JF has to exhibit to provide a viable candidate to reproduce the observed decay behavior of *H*0(*z*) (Equation (36)). We conclude by observing that the picture depicted above relies on the concept of a slow-rolling phase of the scalar field, when it approaches the value *φ* ' 1 and, in this respect, the potential term should have for such value a limiting dynamics, which remains there confined for a sufficiently long phase. It is just in such a limit that we are reproducing a ΛCDM model, but with the additional feature of a slowly varying Einstein constant. As far as the value of *z* increases, the deviation of the considered model from General Relativity becomes more important, but this effect is observed mainly in the gravity-matter coupling. In other words, the motion of the photon, as observed in the gravitational lensing, is not directly affected by the considered deviation, since the geodesic trajectories in the space-time do not directly feel the Einstein constant value. This consideration could allow for a large deviation of *φ* from the unity that is expected in studies of the photons' propagation.
