Signature of f ( R ) gravity via Lemaître-Tolman-Bondi inhomogeneous perturbations

We analyze inhomogeneous cosmological models in the local Universe, described by the Lemaître-Tolman-Bondi (LTB) metric and developed using linear perturbation theory on a homogeneous and isotropic Universe background. Focusing on the diﬀerent evolution of spherical symmetric inhomo-geneities, we compare the Λ LTB model, in which the cosmological constant Λ is included in the LTB formalism, with inhomogeneous cosmological models based on f ( R ) modiﬁed gravity theories viewed in the Jordan frame. We solve the system of ﬁeld equations for both inhomogeneous cosmological models adopting the method of separation of variables: we integrate analytically the radial proﬁles of local perturbations, while their time evolution requires a numerical approach. The main result of the analysis concerns the diﬀerent radial proﬁles of local inhomogeneities due to the presence of a non-minimally coupled scalar ﬁeld in the Jordan frame of f ( R ) gravity. While radial perturbations follow a power-law in the Λ LTB model, Yukawa-like contributions appear in the f ( R ) theory. Interestingly, this latter peculiar behavior of radial proﬁle is not aﬀected by the choice of the f ( R ) functional form. The numerical solution of time-dependent perturbations exhibits a non-diverging proﬁle. This work suggests that investigations about local inhomogeneities in the late Universe may allow us to discriminate if the present cosmic acceleration is caused by a cosmological constant term or a modiﬁed gravity eﬀect.


Introduction
Our current theoretical understanding in cosmology is essentially based on two crucial pillars: 1) General Relativity (GR) is the underlying gravitational theory; 2) the cosmological principle states that the spatial distribution of matter in the Universe is revealed as homogeneous and isotropic on scales sufficiently large.This paper aims to investigate the robustness of these two pillars, examining also possible deviations.
Presently, the cosmological concordance model is generally referred to as the well-known ΛCDM model [1][2][3][4], which involves a cold dark matter (CDM) component and a cosmological constant Λ.This cosmological paradigm is consistent with most of the data, providing a reliable picture of the present-day observed Universe.For instance, maps of the cosmic microwave background radiation (CMB) [5] have corroborated the idea of a homogeneous and isotropic Universe on large scales 1 .The geometry of a purely homogeneous and isotropic Universe is properly described within the framework of the ΛCDM model through the Friedmann-Lemaître-Robertson-Walker (FLRW) metric [3].
Concerning the first pillar of the ΛCDM model, despite all the remarkable predictions and successes to explain theoretically many observational facts in the Universe, it is believed that GR is not the ultimate theory of gravity.Indeed, some open problems in cosmology, such as the need for dark components such as CDM and dark energy, the exact nature of which is still unknown, the cosmological constant problem [7,8], and the Hubble constant tension 2 , may be regarded as a possible signal of the breakdown of GR on galactic and cosmological scales and have motivated people to discuss some modifications of the theory.
Among all the possible proposals of modified gravity with respect to the Einsteinian formulation, the so-called f (R) gravity theories [25][26][27][28][29][30][31][32][33] stand for their simple morphology since de facto only an extra scalar degree of freedom is added to the gravitational field dynamics.These theories represent an extension of GR, a particular class of modified gravity, in which the scalar curvature R in the gravitational Lagrangian density is replaced by a function f (R), i.e., an extra geometrical degree of freedom with respect to GR.Moreover, scalar-tensor gravity provides an equivalent representation of f (R) theories: adopting a formulation in the so-called Jordan frame [26,[28][29][30][31]33], the additional mode is explicitly translated into a non-minimally coupled scalar field to standard gravity 3 .1 This fact is widely accepted, although CMB anomalies and anisotropies have been recently emerged in the form of an observed temperature dipole [6]. 2 In particular, CMB measurements [5] have provided a value of the actual expansion rate, i.e., the Hubble constant, that is incompatible with observations of local probes, such as Cepheids and type Ia supernovae (SNe Ia) [9], with a significance level of 4.9 σ.See [10][11][12][13][14][15] for comprehensive reviews. Se also [16][17][18][19][20][21][22][23][24] for attempts to reduce the tension, using new statistical analysis and combining different cosmological probes.3 As typical in scalar-tensor theories, the f (R) gravity can be also described upon a proper conformal transformation in the so-called Einstein frame [26,[28][29][30]33], which is mathematically equivalent to the Jordan one.Actually, the problem The interest in the f (R) modified gravity is because extra degrees of freedom may allow us to find alternative explanations for the abovementioned unresolved problems in cosmology.For instance, some successful proposals in the f (R) gravity [40][41][42][43] describe the present accelerating Universe, avoiding to introduce ad hoc extra components, such as dark energy.Other examples of f (R) models have been developed to try to alleviate or solve the Hubble constant tension [12,[44][45][46][47][48][49][50].Moreover, the presence of extra degrees of freedom within the f (R) metric formalism is also discussed in gravitational-wave physics, regarding the nature of their polarizations [51].
Therefore, an increasing interest has risen in recent years to develop new methods to discriminate between the standard ΛCDM model and modified gravity cosmological scenarios for the present Universe.Extended f (R) theories of gravity admit a larger number of solutions than Einstein field equations in GR, but extra degrees of freedom and non-linearity imply non-trivial cosmological dynamics.Some cosmological exact solutions in the f (R) gravity have been found for simple scenarios, especially in a homogeneous and isotropic cosmology [29,[52][53][54][55].
Concerning the second pillar of the ΛCDM model, that is the homogeneity and isotropy of the Universe on large scales, it should be emphasized that the cosmological principle is not invariant for any spatial scale: our present Universe seems to reach the conditions to be homogeneous on a scale of about 100 Mpc, as suggested by galaxy surveys and the large-scale structure of the Universe [56,57].Thus, when considering physical phenomena that occur on smaller spatial scales, local features of the Universe lead inevitably to deviations from the FLRW geometry, which could affect cosmological parameters [58][59][60][61][62] and the luminosity distance distribution [63].See also [64][65][66] for the development of the averaging formalism in cosmology to average scalar quantities in a limited region of space-time.
Furthermore, observational evidence of a local void (an underdense region) have emerged on scales of several hundreds of Mpc [67][68][69][70].In this of identifying the physical frame has been long debated [34][35][36][37][38][39].In this work, we decide to adopt the Jordan frame, since it is provided only by a pure redefinition of fields starting from the f (R) metric formalism.
regard, the Lemaître-Tolman-Bondi (LTB) model [4,8,[71][72][73] is widely employed to describe spherically symmetric non-stationary inhomogeneities in the local Universe, while it approaches a homogeneous Universe far enough from the center of symmetry.The LTB solution with the presence of a cosmological constant Λ is commonly referred to as the ΛLTB model, which has been studied for many decades as a possible (simplified) framework to consider local deviations of the Universe today from the homogeneity [74][75][76][77][78][79][80] or to alleviate the Hubble constant tension [81][82][83][84][85].
Since f (R) theories and the ΛLTB model separately can only alleviate the Hubble constant tension, the combination of more than one nonstandard physical effect, i.e., modified gravity and inhomogeneous cosmology, may be needed to accommodate these data inconsistencies into a cosmological model.
However, in this work, we do not focus on the Hubble tension, but we consider a prodromic question, i.e., the formulation of the LTB cosmology in modified gravity.Other papers were carried out to find inhomogeneous LTB solutions in extended gravity scenarios [86][87][88][89][90][91][92]; here we consider the f (R) metric formalism to highlight possible peculiarities of local inhomogeneities compared to the GR solutions.Our aim is to search for some specific markers of the matter distribution, which could allow us to distinguish a standard paradigm from an alternative cosmology.Indeed, local inhomogeneities incorporated in a modified gravitational scenario might become a predictive tool.
In this paper, in the spirit of describing small deviations from the homogeneity of the Universe today, we treat the inhomogeneous LTB metric as the flat FLRW background solution with the addition of small spherically symmetric perturbations, following a linear perturbation approach 4 .We extend the work carried out in [96,97].The present analysis aims to discriminate between the ΛLTB cosmological paradigm and the LTB solution as emerging in the Jordan frame of the f (R) gravity, comparing the different evolution of inhomogeneous perturbations within these two schemes to describe the local behavior of the actual Universe.In this regard, we evaluate the role of a cosmological constant with respect to the presence of a non-minimally coupled scalar field in the Jordan frame, examining separately the cosmological dynamics.In particular, regarding the background cosmology in the Jordan frame, we adopt the f (R) model developed by Hu and Sawicki [40], which provides an interesting alternative to the dark energy component for the current cosmic acceleration within the f (R) gravity.
Therefore, we investigate the first-order perturbation field equations for local inhomogeneities both in the ΛLTB and f (R) inhomogeneous cosmological models.More specifically, the analysis of the 0-1 component of the gravitational field equations points out discrepancies between the two formalisms considered, due to the presence of the non-minimally coupled scalar field.In our analysis of the first-order perturbation equations, we adopt the separable variables method as a mathematical technique to address the solution of the partial differential equations system.Then, we separately obtain the time evolutions and radial profiles of perturbations within both two cosmological models.Actually, we get an analytic expression only for the radial part, while we need to numerically evaluate the time evolution of perturbations.
The main result of our work is the different radial patterns of inhomogeneous perturbations in GR and the f (R) cosmology.Indeed, in the former case, we deal with a radial solution that is basically a power law; in the latter case, a peculiar Yukawa-like dependence emerges in the Jordan frame gravity whatever the f (R) model considered.Both these radial solutions suggest that inhomogeneities decay rapidly on large scales according to the cosmological principle, but furthermore, from a theoretical point of view, we have a specific trace of how a modified gravity model can be distinguished from a standard gravity scenario in the presence of a cosmological constant.In other words, the Yukawa-like radial solution is a peculiar feature of the f (R) gravity.
Regarding our numerical analysis of the time dependence of the obtained solutions, we show that local inhomogeneities are governed by a dynamics preserving the stability of the isotropic Universe both in the ΛLTB and modified gravity theories.Therefore, these stable-time solutions are physically acceptable and can be thought of as large-scale corrections to the FLRW geometry.
This work provides an interesting arena to investigate different scenarios of the clumpy inhomogeneous cosmology; at the same time, the method developed in this paper supplies a useful criterion to distinguish between GR and f (R) modified gravity models on cosmological scales.
This paper is structured as follows: in Sect. 2 we introduce the f (R) modified theories of gravity in the Jordan frame and the Hu-Sawicki model; in Sect. 3 and Sect. 4 we implement the cosmological dynamics to the ΛLTB model in GR and the LTB spherically symmetric solution as emerging in the Jordan frame of the f (R) gravity, respectively; in Sect. 5 we show our perturbation approach to study local inhomogeneities in the ΛLTB model, obtaining background and first-order perturbation solutions, while in Sect.6 we proceed similarly for the f (R) theories in the Jordan frame; lastly, we summarize our results and conclusions in Sect.7.
We adopt the metric signature (−, +, +, +) throughout the paper, and we use natural units for the speed of light c = 1.We denote with χ ≡ 8 π G the Einstein constant, being G the gravitational Newton constant.The components of the gravitational field equations are sometimes denoted with a pair of indices µ − ν, in which the indices µ, ν = 0, 1, 2, 3 move along the set of coordinates.

f (R) modified gravity in the Jordan frame
We investigate the equivalence between f (R) theories of gravity and the scalar-tensor representation in the Jordan frame [25][26][27][28][29][30][31][32][33].There is no a priori motivation to consider a linear gravitational Lagrangian density with respect to the Ricci scalar R, apart from obtaining a second-order partial differential system for the field equations.The gravitational Lagrangian density, which is given by R in GR in the Einstein-Hilbert action, is generalized in the context of f (R) theories as a function f of R, i.e., an extra degree of freedom.In addition, considering a matter term S M , the total action of f (R) gravity is given by where g is the determinant of the metric tensor with components g µν , and ψ refers to the matter fields.
It can be shown that the extended gravitational field equations within the f (R) metric formalism are fourth-order partial differential equations in the metric.If f (R) = R, specifically, the fourth-order terms vanish, and field equations reproduce exactly the Einstein field equations in GR.
To bring field equations to a form that is easier to handle, the f (R) gravity can be restated in the scalar-tensor formalism.In particular, in the socalled Jordan frame 5 , it can be checked that the following action is dynamically equivalent to the f (R) action given by Eq. ( 1) if f ′′ (R) = 0, where the scalar field is governed by the scalar field potential defined as It should be noted that the extra degree of freedom given by the f (R) function turns into a scalar field 6 φ, which is non-minimally coupled to the metric.Conversely, there is only a minimal coupling for the matter.Variations of the action (2) with respect to the metric and scalar field lead to field equations in 5 It is quite similar to the prototype of extended gravitational theories, the Brans-Dicke formulation [98,99] with a non-zero potential and a null Brans-Dicke parameter or the O'Hanlon proposal [100]. 6Actually, as shown in [101], the relation f ′′ (R) = 0 is a redundant requirement for the dynamic equivalence between an f (R) theory and a scalar-tensor formulation in the Jordan frame.Indeed, it is sufficient to require that f ′ (R) be invertible (continuous and one-to-one in a given interval), i.e., the existence of R = R f ′ .In this way, it is possible to build a scalar field potential V (φ).
the Jordan frame respectively, where G µν is the Einstein tensor, T µν is the stress-energy tensor of matter, ∇ µ is the covariant derivative associated with the Levi-Civita connection of the metric, and Furthermore, by taking the trace of Eq. (5a) and using Eq.(5b), a dynamical equation for the scalar field is obtained for a given matter source, where T is the trace of the matter stress-energy tensor.
Although the presence of a non-minimally coupled scalar field implies non-trivial dynamics in the Jordan frame, it is often convenient to adopt the Jordan frame of the f (R) gravity, since the field equations (5a), (5b), and ( 6) are now secondorder differential equations.
Looking at the action given in Eq. ( 2), the extra degree of freedom provided by f (R) does not affect the matter action even in the Jordan frame.Therefore, the stress-energy tensor of ordinary matter must be divergence-free as in GR: ∇ ν T µν = 0. Since Bianchi identities are kept in modified gravity, the Einstein tensor is also covariant divergence-free: ∇ ν G µν = 0. Here, we define an effective stress-energy tensor related to the scalar field (7) to rewrite the field equations (5a) in the Jordan frame as Finally, considering the divergence-free relations above, it is trivial to show that The effective stress-energy tensor T [φ] µν does not satisfy the usual law of the ordinary matter [102][103][104][105], unless in vacuum (T µν = 0).The extra scalar degree of freedom in the Jordan frame is quite different from a matter field; actually, φ is an effective scalar field originating from a scalar mode intrinsically due to a modification in the gravitational action.Note that the laws (9) describing the dynamics of T µν , as well as the continuity equation related to T µν , are not independent of the field Eqs.(5a) and (6).However, Eqs. ( 9) can be employed as auxiliary equations to rewrite field equations in a different equivalent form.

The Hu-Sawicki model
Among several proposals for the functional form of the f (R), one of the most studied dark energy models is provided by Hu and Sawicki [40,41].It is useful to refer to the deviation F (R) from the gravitational Lagrangian density in GR, i.e.
The functional form of the deviation F (R) for the Hu-Sawicki (HS) model is where n is a positive integer, c 1 and c 2 are the HS dimensionless parameters, and m 2 ≡ χ ρ m0 /3 with ρ m0 matter density today.
Note that for R ≫ m 2 the limiting case with an effective cosmological constant Λ eff = c 1 m 2 /2 c 2 is recovered.Approximating the cosmic accelerated phase of a flat ΛCDM model with an effective cosmological constant, we obtain the first constraint on the parameters c 1 and c 2 : in which we have used the definitions of cosmological density parameters Ω m0 ≡ ρ m0 /ρ c0 and Ω Λ0 ≡ ρ Λ0 /ρ c0 for the matter component and the cosmological constant, respectively, being ρ c0 ≡ 3H 2 0 /χ the critical energy density of the Universe today, ρ Λ0 ≡ Λ/χ the energy density associated with Λ, and H 0 is the Hubble constant.The subscript 0 denotes the present cosmic time t 0 at the redshift z = 0. Furthermore, Hu and Sawicki [40] have shown that today R 0 ≫ m 2 , and also the approximation R ≫ m 2 is viable for the entire past cosmic expansion.
Hereinafter, we set n = 1 for simplicity, focusing on the most extensively studied scenario in the HS gravity.In that case, the derivative In particular, for a flat ΛCDM model where a = a (t) is the scale factor in terms of the cosmic time.We have used the well-known Friedmann equations and the definition of the Ricci scalar R in a flat FLRW metric [3].We can rewrite Eq. ( 13) evaluated today in the limiting case for R ≫ m 2 : in which we have implicitly assumed that the value of Λ eff to be the same as Λ in the ΛCDM limit.
We have also set the conventional notation for the scale factor today a 0 = 1.Hence, by setting a reference value for F R0 , we can obtain the second constraint on c 1 and c 2 .Note that, according to Eqs. ( 3) and ( 10), we have in the Jordan frame: It should be noted that F R quantifies the deviation from the GR scenario, where φ = 1.
The scalar field potential V (φ) in the Jordan frame for the HS model assumes the following form: where we used the definitions given in Eqs. ( 3) and ( 4) referred to the F (R) function in Eq. (11).Moreover, we have selected the branch for the potential related to a minus sign just before the Scalar field ϕ Hu-Sawicki V(ϕ)/m 2 Fig. 1 Behavior of the HS scalar field potential in the Jordan frame, defined in Eq. ( 17).It should be noted that square root in Eq. ( 17) to converge to an asymptotically stable de Sitter Universe [106,107].
Finally, we show two reasonable values for the HS dimensionless parameters c 1 and c 2 .More precisely, we fix Ω m0 = 0.3111, and Ω Λ0 = 0.6889 from the Planck measurements [5]; we set the value of the derivative of the field at the present cosmic time |F R0 | = 1.0 × 10 −7 , considering the strongest bound between solar system [40] and cosmological constraints [108][109][110].Then, using the conditions ( 12) and (15) for the ΛCDM limit, we obtain c 1 = 2.0 × 10 6 and c 2 = 1.5 × 10 5 .Considering these values for c 1 and c 2 , in Fig. 1 we show the profile of the HS scalar field potential in the Jordan frame, noting a slow evolution of V (φ).
In the LTB model, the space is isotropic only observing the Universe from a specific preferred point, i.e., the center that is singled out by adopting a spherical symmetry, where an observer is supposed to be located.Hereinafter, we consider the evolution of a dust cosmological model resulting in a spherical mass overdensity (or underdensity) with vanishing pressure (p = 0), which can be formulated in the LTB formalism.Furthermore, we focus on the late Universe, hence we neglect relativistic species, since they are subdominant today, i.e., Ω r0 ∼ 10 −5 .
The LTB spherically symmetric line element in the synchronous gauge is written as in which t is the cosmic time, r is the radial coordinate indicating the spatial distance from the preferred point, and dΩ is the solid angle element.Note that there are two metric functions: α = α (t, r) and β = β (t, r).
The three independent Einstein field equations in the ΛLTB model with a pressure-less dust are provided by the 0-1, 0-0, and 1-1 components, which rewrite as respectively, where () = d/dt and () ′ = d/dr.We have separated explicitly the cosmological constant and matter term, and ρ is simply the energy density of the matter component (we neglected the subscript m in ρ for brevity).The other non-null field equations in the LTB metric are related to the previous equations system due to the spherical symmetry.More precisely, it is straightforward to show the following relations between the Einstein tensor components in the LTB metric: and also It should be noted that the two metric functions α and β can be related in GR, by exploiting the 0-1 component (19a) of the Einstein field equations, in order to rewrite the LTB line element (18) in a simpler form [111]. Indeed, Eq. (19a) admits the solution where g (r) is an arbitrary function of the radial coordinate r.Then, the LTB line element (18) can be rewritten as in which the following parametrization has been adopted with K = K (r), and a (t, r) ≡ e β r −1 has been defined as the generalization of the scale factor in an inhomogeneous Universe.By using the LTB metric, the remaining Einstein field equations (19b) and (19c) become respectively.It may be observed that the form of the LTB metric given by Eq. ( 22) reminds the FLRW line element.More specifically, if a (t, r) and the LTB curvature function K (r) do not depend on the radial coordinate r, the FLRW geometry is exactly recovered to describe a homogeneous and isotropic Universe.Furthermore, it is straightforward to show in that limit that Eqs.(24a) and (24b) turn into the Friedmann equations in the FLRW metric.
Finally, the continuity equation for a pressureless perfect fluid in the LTB metric (18) can be written as or equivalently if the LTB metric in the form given by Eq. ( 22) is considered.We recall that the energy conservation law results from the divergenceless law of the stress-energy tensor ∇ µ T µν = 0 for ν = 0.
4 The LTB model in the Jordan frame of f (R) gravity We study the cosmological dynamics in the LTB metric (18) within the framework of the f (R) gravity in the Jordan frame.The 0-1, 0-0, 1-1 components of the gravitational field equations (5a) are written as β′ respectively, where we have put a cosmological pressure-less dust as a source.In the Appendix A, it is shown that the other non-vanishing gravitational field equations, i.e., 2-2 and 3-3 components, depend on the previous set of equations, as it must be, basically due to the spherical symmetry in the LTB geometry.Moreover, the scalar field equation ( 6) rewrites as Note that in an inhomogeneous cosmology all the quantities φ, ρ, α, and β depend on both t and r.
It should be emphasized the occurrence of extra contributions in field equations ( 27) and (28) with respect to the ΛLTB model, due to the coupling between the scalar field φ and the metric functions α, β, as well as the presence of the scalar field potential.For instance, it is quite clear to recognize an extra coupling term by comparing the 0-1 field Eq. (27a) with the respective Eq. (19a) in GR.As a consequence, this coupling in the Jordan frame does not allow us to find a relation between α and β, and then rewrite the LTB metric in a simpler form, unlike the ΛLTB model in GR (Sect.3).For all these reasons, the cosmological dynamics in the Jordan frame is really different from the GR scenario.
On the opposite, the continuity equation related to the ordinary stress-energy tensor T µν for a dust in the LTB metric exactly exhibits the same form provided in Eq. ( 25) both in GR and f (R) gravity.In the latter theory, other additional equations are those related to the effective stress-energy tensor T [φ] µν , defined in Eq. ( 7), for the scalar field in the Jordan frame.More specifically, Eqs.(9) in the LTB metric become It should be recalled that these laws for T [φ] µν are not independent of the field equations ( 27) and (28), because these laws basically come from field equations using Bianchi identities.Nevertheless, these additional equations ( 29) can be useful to rewrite field equations in a different form.For instance, Eq. (29b) has been employed to find the dependence between the 1-1 and 2-2 components of the gravitational field equations (27c) and (A2), which implies spatial isotropy in the LTB geometry, as it has been shown in the Appendix A.
To sum up, within the Jordan frame of f (R) gravity in the LTB metric, we have obtained a system of four partial differential equations ( 27) and ( 28) with four unknown functions: α (t, r), β (t, r), ρ (t, r), and φ (t, r).Note that the scalar field potential V (φ) provides a degree of freedom in the theory.Furthermore, other supplementary equations are provided by Eqs. ( 25) and (29).

Perturbation approach for the LTB model in General Relativity
In this section, we consider local inhomogeneities of the Universe as small spherically symmetric perturbations over a flat background FLRW geometry.We follow a linear perturbation approach, so that we have the FLRW geometry at the zerothorder perturbation theory, while we build a lumpy Universe described by the LTB metric at the firstorder perturbation [96].Thus, we can write the LTB metric tensor components as Hereinafter, we use an overbar to denote quantities referred to the background homogeneous and isotropic Universe, and the symbol δ is related to linear perturbation terms.We emphasize that, choosing this decomposition in Eq. ( 30), we require spherically symmetric perturbations.Moreover, we adopt the synchronous gauge for the LTB metric, which intrinsically includes two degrees of freedom in the perturbed metric.Indeed, two independent metric functions, i.e., α (t, r) and β (t, r), are contained in the original LTB line element (18) or, equivalently, a (t, r) and K 2 (r) in GR, according to Eq. ( 22).Actually, K 2 (r) is not exactly a dynamical degree of freedom but an arbitrary parametric function, as a result of the LTB cosmological dynamics.To be more specific, we recall that, in Sect.3, the 0-1 component (19a) of the Einstein field equations has allowed us to find a relation between the metric functions α (t, r) and β (t, r) in GR, hence to reduce one degree of freedom, and the LTB metric in GR assumes the form given in Eq. ( 22) in terms of a (t, r) and K 2 (r).
Since the background FLRW and the perturbed LTB metrics are both locally rotationally symmetric and are given in the same normal geodesic frame, we only need to focus on scalar functions, as shown in [112,113].Hence, the scale factor a (t, r) and the energy density of the matter component ρ (t, r), contained in the field equations (24a) and (24b), are defined as: Note that the curvature function K 2 (r) in the LTB metric (22) involves radial inhomogeneities, so it can be regarded as a perturbative contribution.In other words, even if hereinafter we set the curvature parameter in the FLRW metric as k = 0, it is possible to define a linear curvature pertubation, which is exactly given by the curvature function K 2 (r) in the LTB metric, according to the metric decomposition in Eq. (30).Thus, in addition to the metric perturbation δa, we have another perturbed quantity with respect to the background FLRW metric, which is K 2 (r).
We also require that background terms dominate over the linear perturbations at any time t, or redshift z, in the late Universe.This condition is dictated by the cosmological principle.
Once the decomposition has been defined in Eqs.(31a), (31b), the evolution in time and space of the physical quantities can be obtained by studying the gravitational field equations (24a) and (24b) at background and linear levels.
In particular, we define a dimensionless time variable in which t 0 is the present cosmic time (today τ = 1).Note that τ is the cosmic time in units of the present Hubble time since t 0 can be approximately written in terms of the Hubble constant as t 0 ≈ 1/H 0 [3].

Background solution
If we consider only background terms, it is straightforward to show that Eq. (24a) turns into the first Friedmann equation in a flat FLRW geometry: being H(t) the Hubble parameter.Furthermore, the other Eq.(24b) can be rewritten at the background level, and combined with Eq. ( 34), provides the second Friedmann equation, also commonly named the cosmic acceleration equation: Then, we rewrite the first Friedmann equation (34) in terms of τ , defined in Eq. ( 33), as We The Friedmann equation ( 36) admits an analytical solution in the late Universe, that is the background scale factor in terms of τ : Note also from Eq. ( 37) that the deceleration parameter q (τ ) ≡ − ä ā−1 H −2 → −1 for τ → +∞, as it should be in a dark-energy-fullydominated Universe.We here stress that the solution (37) applies only in the late Universe, otherwise we need to solve numerically the field equations, if we also consider the radiation contribution.Nevertheless, we are interested in the evolution of local inhomogeneities at late times.
It should be noted that in the limit we recover the well-known relation for the background scale factor ā satisfied in the matterdominated Universe [3].
Concerning the evolution of the background energy density ρ, we start from the continuity equation (26) in the LTB metric.We verify that this equation at the zeroth order becomes simply the respective continuity equation for the matter component in the FLRW metric, i.e.
which provides ρ ∼ ā−3 .Actually, we recall that Eq. ( 39) is not independent of the two Friedmann equations ( 34) and (35), since it can be derived by combining them.Finally, we focus on the relation between the redshift z and the dimensionless parameter τ , defined in Eq. (33), to understand to what extent of τ values the background solution (37) can be applied in the late Universe.To be more accurate in this computation, we also include the radiation contribution to write τ as where we have used the Friedmann equation (36), the definition of τ in Eq. ( 33), and we recall that ā0 /ā = 1 + z.The quantity τ (z) can be computed numerically for a given redshift z, after specifying the cosmological parameters: Ω m0 = 0.3111, Ω Λ0 = 0.6889, and Ω r0 = 9.138 × 10 −5 from Table 2 in [5].For instance, we can compute: τ (z eq ) = 0.046 at the redshift of the matterradiation equality z eq = 3403.5;τ (z DE ) = 0.75 at the matter-dark energy equality z DE = 0.303; τ (z 100eq ) = 0.052 when the energy density of the matter component was one hundred times more than the radiation contribution at z 100eq = 33.045.In particular, considering this latter value of τ , you can see in Fig. 2 the behavior of the background scale factor ā (τ ) and the deceleration parameter q (τ ) in the range τ 100eq < τ < 5, when relativistic species are negligible.

Linearly perturbed solutions in an inhomogeneous Universe
We analyze the impact of local inhomogeneities in the cosmological dynamics.The linearized field equations allow us to investigate the evolution of spherically symmetric perturbations.
If we include local inhomogeneities in the firstorder perturbation theory, the Eqs.Fig. 2 Evolution of the background scale factor ā (τ ) (top panel) and the respective deceleration parameter q (τ ) (bottom panel) in terms of the time dimensionless parameter τ , defined in Eq. ( 33), for a flat ΛCDM model within the range τ 100eq < τ < 5, according to the solution given in Eq. (37).
respectively.Furthermore, we rewrite the continuity equation ( 26) at linear order: To study separately the evolution of local inhomogeneities in time and space, we adopt the separation of variables method to solve analytically the first-order perturbation equations.Hence, we define time and radial functions for all linear perturbations: (43) The quantities а p (t) and a p (r) are both dimensionless.We assume, without loss of generality, that R p (t) has the physical dimensions of an energy density as ρ and δρ, while we treat ̺ p (r) like a dimensionless quantity.We also recall that the curvature perturbation K 2 (r) in the LTB metric depends only on the radial coordinate.
We stress that if we would also include nonlinear terms, then the separation of variables could not lead to a general solution.However, the linearization procedure adopted for the dynamics allows us to use a separation of variables characterized by the factorization (43) of the time and space dependences in the linear perturbation theory.
Using the factorization (43), Eq. (41b) can be split into two parts.By setting the radial dependence as we obtain an ordinary differential equation for the time evolution: in which we have also considered that ā ≫ а p .We would emphasize that the assumption given in Eq. ( 44) is suggested by the form of Eq. ( 41b), once we used the separation of variables from Eq. ( 43).Nevertheless, we have still two metric perturbations given by the quantities а p (t) and a p (r).
We proceed similarly for the first-order perturbation continuity equation (42).After straightforward calculations, by using again the factorization (43) in Eq. ( 42), we separate terms that depend only on t from those related to r.In particular, the time evolution is provided by while we obtain the following radial dependence The constant X is introduced by the separation of variables method.We recall that the background solution for ā is written in Eq. ( 37) and also ρ ∼ ā−3 .We rewrite the first-order perturbation Eq. (41a), by employing Eq. ( 43) and dividing both sides of the equation by ā2 а p a p , as We observe the presence of several mixed terms depending both on t and r.However, we can reduce the number of these mixed terms, by using Eqs.( 34), (44), and (47) to rewrite Eq. ( 48) as Since we want to avoid trivial solutions, we should have ̺ p = 0 and a p = 0.Then, we obtain a single ordinary differential equation in the time domain: It is straightforward to check the compatibility between Eqs. ( 45), (46), and (50).Indeed, by combining the time derivative of Eq. ( 50) with Eq. ( 46), the background field equations (34), and (35), it is easy to build exactly Eq. (45).
Then, we rewrite the term in the brackets on the right-hand side of Eq. ( 46) by using Eq. ( 50), and we get an ordinary differential equation with a single variable R p (t), that is Therefore, focusing on the time domain, we can obtain numerical solutions for the unknown quantities а p and R p from the linearized equations ( 45) and (51).In particular, Eq. ( 45) shows exactly the same behavior in terms of t and τ , as it can be checked by using the definition of τ (33).Recalling the expression (37) of the background scale factor in GR, Eq. ( 45) can be solved numerically.We set the initial conditions at τ = 1 today: а p (1) = 10 −5 and аp (1) = 0.Moreover, we fixed the same values for Ω m0 and Ω Λ0 adopted in Sect.5.1.In the upper panel of Fig. 3, you can see the numerical results for 1 ≤ τ ≤ 5. Note that the perturbed scale factor а p increases as τ grows, and this fact may be a problem if perturbations become unstable.However, the evolution of а p is dominated by the background term ā at any time τ .Indeed, the ratio between the perturbation and background terms with η (τ ) ≡ |а p /ā| ≪ 1 for any τ , as it is shown in the middle panel of Fig. 3.In other words, perturbations of the scale factor due to local inhomogeneities will remain small over time.
Concerning the time evolution of the perturbed energy density of the matter component R p , we rewrite Eq. ( 51) in terms of τ as in which we have defined the dimensionless quantity O p ≡ R p / ρc0 .We solve numerically Eq. ( 52) with the initial condition O p (1) = 10 −5 for τ = 1.
The results are shown in the bottom panel of Fig. 3.Note that O p (τ ) ≪ 1 for 1 ≤ τ ≤ 5. Finally, after obtaining numerical solutions for the time domain, we focus on the radial part of perturbations.We recall that K 2 (r) = a p (r) from Eq. ( 44).Moreover, we assume for simplicity a proportionality between radial perturbations, i.e., ̺ p = C a p with a constant C.
Then, we solve the ordinary differential equation ( 47) in r, and we obtain a power-law behavior for the radial correction of the perturbed scale factor in GR: for which y ≡ 3 − C X.More in detail, we need to impose the condition y > 0 to ensure that inhomogeneities decay on large scales, according to the cosmological principle.
6 Perturbation approach for the LTB model in the Jordan Frame of f (R) gravity In this section, we compare and discuss the evolution of inhomogeneous perturbations in GR and in the Jordan frame of f (R) gravity, considering again a cosmological dust (p = 0) in the LTB geometry.A complete general solution in linear perturbation theory within the metric f (R) gravity was developed in [40,41,114].Instead of proceeding with a fourth-order cosmological ) Fig. 3 Top panel: evolution of the linearly perturbed scale factor аp (τ ) in units of 10 −4 and in terms of the parameter τ , provided by the numerical solution of Eq. ( 45).Middle panel: the ratio between the first-order perturbation term and background scale factor η (τ ) ≡ |аp/ā| versus τ in units of 10 −5 .Bottom panel: evolution of the dimensionless perturbed energy density Op (τ ) ≡ Rp/ρ c0 in units of 10 −6 .Note that all perturbed contributions are smaller than respective background terms for 1 ≤ τ ≤ 5.
dynamics, here we work in the equivalent Jordan frame f (R) gravity.Furthermore, as a particular case, we focus on spherically symmetric perturbations, following the same perturbation approach developed in Sect. 5 by using the metric decomposition in Eq. ( 30).Hence, we split the metric functions α and β, the energy density ρ, and the scalar field φ into background terms plus linear corrections as We also require again that inhomogeneities are much smaller than respective background terms.
It should be stressed that, in the Jordan frame, we can not use the LTB metric in the simpler form (22), but we refer to the original LTB line element (18).Thus, the two degrees of freedom of the perturbed metric are given by α (t, r) and β (t, r); we no longer refer to a (t, r) and K 2 (r), as in the GR scenario.
Note that the background quantities ᾱ and β are related to the scale factor ā, since we want to reproduce a flat FLRW geometry at the zerothorder perturbation.Then, comparing a flat FLRW metric with the LTB line element in the form given by Eq. ( 18), it is straightforward to show that ᾱ (t) = ln (ā (t)) , β (t, r) = ln (ā (t) r) .
(55) As a consequence, the metric tensor component g rr in the LTB metric can be approximated for δα ≪ 1/2 as g rr = e 2α = e 2( ᾱ+δα) ≈ ā2 (t) (1 + 2 δα) , (56) in which it is possible to recognize the respective metric tensor component g rr in a flat FLRW metric as background term.We can follow the same reasoning for the other metric tensor components involving β with the assumption δβ ≪ 1/2.Also, we need to expand the scalar field potential V (φ), defined in Eq. ( 4), which appears in the field equations (27b), (27c), (28).Therefore, including local inhomogeneities, we obtain for the first-order perturbation theory.Similarly, we can rewrite the derivative of V (φ).
As we have proceeded to study the gravitational field equations in GR in Sect. 5 to find background and linear solutions, now we analyze the dynamics of f (R) gravity in the Jordan frame, provided by the field equations (27), and (28) in the LTB metric.

Background solution
If we consider all background quantities, we do not include spherically symmetric perturbations and Eqs.(27b), (27c), and (28) turn into field equations in the Jordan frame of f (R) gravity in a flat FLRW geometry: Eq. ( 58a) is the modified Friedmann equation, Eq. ( 58b) is the modified acceleration equation, and Eq.(58c) concerns the scalar field evolution in the Jordan frame.In particular, Eq. ( 58b) is obtained by combining the zeroth-order perturbation Eqs. ( 27b) and (27c).Furthermore, Eq. (27a) vanishes at the background level in the synchronous gauge, and it becomes a trivial identity.
Then, we focus on the f (R) HS model in the Jordan frame, which has been introduced in Sect.2.1.Recalling the form of the scalar field potential V φ in Eq. ( 17), we observe that the field equations (58) do not admit any analytical solutions, and we have to solve it numerically.
In this regard, we rewrite the full set of equations (58) in terms of the dimensionless parameter τ , defined in Eq. ( 33), and we obtain, respectively: ) Fig. 4 Numerical background solutions, considering a flat FLRW geometry, for the f (R) HS model in the Jordan frame.Top panel: evolution of the scale factor ā in terms of the parameter τ .Middle panel: behavior of the deceleration parameter q (τ ).Bottom panel: the deviation from the GR scenario ( φ = 1) with units of 10 −7 for the vertical axis.
We choose the parameters of the model in such a way that the background modified gravity scenario is almost equivalent to the ΛCDM cosmological model with the aim of focusing later on the differences between the linear perturbation solutions in the two models.More precisely, we fix Ω m0 = 0.3111, the same value adopted in Sect.5.1, and we set the value |F R0 | = 1.0 × 10 −7 at the present cosmic time (redshift z = 0 or τ = 1), which provides information about the deviation from the GR scenario, according to Eq. ( 16).As a consequence, we constrain the HS dimensionless parameters: c 1 = 2.0 × 10 6 and c 2 = 1.5 × 10 5 , as developed in Sect.2.1.We recall that the profile of the background quantity V φ /m 2 is plotted in Fig. 1.
To guarantee a nearly frozen evolution of the scalar field φ for increasing τ , we impose the following condition: d φ dτ (τ = 1) = 0. Finally, we solve numerically Eqs.(59a) and (59c) for 0.2 < τ < 5, when the relativistic components remain negligible as compared with the matter, to obtain the evolution of ā (τ ) and φ (τ ).The numerical results are shown in Fig. 4. By comparing it with Fig. 2, it should be emphasized that ā (τ ) and q (τ ) exhibit almost the same behavior in GR and in the Jordan frame of f (R) gravity, as desired according to the choice of model parameters abovementioned.In particular, in the bottom panel of Fig. 4, we plot the quantity 1 − φ (τ ) to evaluate the deviation from GR ( φ = 1): we observe more relevant deviations in the late Universe for τ > 1, but nevertheless the background modified gravity dynamics still remains almost undistinguishable from the ΛCDM scenario.Hence, we can shift the attention towards the first-order perturbation solutions in the Jordan frame of f (R) gravity.

Linearly perturbed solutions in an inhomogeneous Universe
We focus on the first-order perturbed equations to study the evolution of spherically symmetric perturbations.Considering the split between background terms and linear perturbations according to Eq. ( 54), the set of field equations (27) becomes Similarly, starting from Eq. ( 28), the linearized scalar field equation is given by Moreover, the continuity equation ( 25) rewrites as at the linear perturbation order.We have used Eq. ( 55) to rewrite ᾱ and β in terms of ā.
Following the same approach we adopted in GR in Sect.5, we use a separation of variables method to study separately the evolution of inhomogeneities in time and space.Hence, we factorize the linear perturbations as: Assuming this factorization, we can rewrite Eqs. ( 60), (61), and (62).However, we notice the presence of several mixed terms depending both on t and r, which do not allow us to solve the equations using the separation of variables in a standard way, unless we rely on reasonable and simplifying assumptions (further details on explicit calculations are in the Appendix B).For instance, we are able to use the separation of variables method for all field equations, if we require the two following conditions: where λ 1 and λ 2 are two proportionality constants.These conditions allow us to simplify the equation system and easily separate time and radial dependences: we obtain a set of differential equations describing the radial profiles of perturbations and another equation system concerning only the time evolution.Hence, starting from Eqs. ( 60), (61), and ( 62), after long but straightforward calculations (see the Appendix B), we obtain a set of equations for the radial part: where µ 1 , µ 2 , µ 3 , and µ 4 are constants arising from the separation of variables.We have four unknown quantities (A p , B p , ̺ p , and ϕ p ), which are related through the condition given by Eq. (64b), the perturbed 0-1 component of field equations (65a), the perturbed 1-1 component (65b), the linearized scalar field equation (65c), and in addition the perturbed continuity equation (65d).In particular, note that µ 3 has dimensions of reciprocal length, i.e., [µ 3 ] = L −1 , as you can see from Eq. (65c).Similarly, in the Appendix B, we write an equation system for the time evolution of perturbations: We have four unknown quantities (A p , B p , P p , and Φ p ) for the time evolution of inhomogeneities, which are fully described by the assumption (64a), the perturbed 0-1, 0-0, 1-1 components of field equations given by Eqs.(66a), (66b), (66c), respectively, the linearized scalar field equation (66d), and the perturbed continuity equation (66e).Note that two equations of the latter list are redundant since the scalar field and continuity equations are not independent of the other field equations.
It should be stressed that the potential V (φ) affects only the time evolution.As a consequence, since the scalar field potential is related to a specific modified f (R) model, the time evolution strongly depends on the particular modified gravity model considered, while the radial part of the linearized equations is completely model free.This is a crucial point to identify a peculiar feature of the inhomogeneities evolution in the Jordan frame of f (R) gravity through analysis of the radial profiles of perturbations.
Once we have split all field equations in time and radial contributions, we seek linear order solutions separately.

Radial profiles
If we focus on the radial evolution of perturbations, we can solve the respective equation system analytically.More specifically, in solving the differential equation (65c), we obtain the Yukawa behavior for the radial solution of the linearly perturbed scalar field: We have set perturbations to vanish at infinity according to the cosmological principle, hence we have only one integration constant γ.We recall that [µ 3 ] = L −1 , as it can be checked also from Eq. (67).Note that ̺ p (r) has the same radial dependence, i.e.
in which we considered the assumption given by Eq. (64b).Consequently, from Eq. (65b), we also get where we used the solution (67).Finally, we combine Eqs.(65a), (67), and ( 69) to obtain the scalar perturbation (70) It should be noted that, in order to satisfy the compatibility between Eqs. (65a), (64b), (65c), (65b), and (65d), we require the following relation between the constants involved in the set of equations 7 .In particular, note that [µ 4 ] = L −2 , since µ 3 has the dimension of inverse length, and the other constants are dimensionless.
It should be pointed out that the f (R) gravity establishes a typical radial scale, r c ≡ µ −1 3 , such that local inhomogeneities vanish for r ≫ r c more rapidly than the respective perturbations in the ΛLTB model.For instance, this fact can be shown by comparing the behavior of the perturbed scale factor given in Eq. ( 53) with the radial profiles of the LTB metric functions in Eqs. ( 69) and (70).It is important to stress that our result is completely independent of the choice of scalar field potential V φ .Hence these radial solutions apply to any f (R) extended model, providing a remarkable feature of the radial evolution within the Jordan frame of the f (R) gravity as compared to GR.

Time evolution
Now we focus on the time evolution of perturbations in the Jordan frame.Clearly, the choice of the background f (R) modified gravity model affects the time dependence of inhomogeneities since V φ appears in almost all equations regarding the time part.Our main result concerns the peculiarity of the radial profiles of perturbations in the Jordan frame, which does not depend on a specific modified gravity model.Here, we merely want to prove the existence of at least one stable time solution.
In particular, Eqs.(64a) and (66a) maintain the same form if we write them in terms of τ as respectively.Furthermore, the linearized scalar field equation (66d) becomes We have defined the density contrast Ξ (τ ) ≡ P p (τ ) / ρ (τ ) for linear perturbations.Note also that the quantity μ3 ≡ µ 3 /H 0 is dimensionless in natural units.Similarly, Eq. (66b) in terms of τ rewrites as and Eq.(66c) turns into: In the latter equation, we have introduced the dimensionless quantity μ4 ≡ µ 4 /H 2 0 in natural units.Finally, considering the definition of the density contrast Ξ (τ ) and the background equation (39), it is straightforward to show that the linearly perturbed continuity equation (66e) simply rewrites in terms of τ as Note that we can use Eq. ( 77) to reduce the number of unknown quantities in the other field equations ( 74).Now we want to solve numerically Eqs. ( 73) and ( 74) to obtain Ξ (τ ) and Φ p (τ ).In this regard, we fix the values of constants for simplicity µ 1 = 1, µ 2 = 10 2 , λ 2 = 5, μ3 = 1 and set initial conditions at τ = 1 for inhomogeneities: Ξ (1) = 10 −5 , Φ p (1) = 10 −12 , and dΦ p /dτ (1) = 0.These values have been chosen to have weak and slowly varying perturbations and guarantee the existence of a stable numerical solution.Moreover, we use the background numerical solutions for ā (τ ) and φ (τ ), and the scalar field potential V φ for the HS model discussed in Sect.6.1.
Then, we obtain numerical solutions for Ξ (τ ) and Φ p (τ ), which are plotted in Fig. 5.It should be noted that these solutions are stable in time, since inhomogeneous perturbations are dominated by background terms for any τ .
Once we have obtained Ξ and Φ p numerically, the remaining perturbed quantities A p and B p can be easily found by using Eqs.( 72) and (77).

Conclusions
In this paper, we have analyzed the LTB spherically symmetric solution, linearized over a flat FLRW background, comparing its morphology in GR and f (R) modified gravity theories, as viewed in the Jordan frame.In the former case, we have referred to the ΛLTB model, including a matter fluid and a cosmological constant; in the latter model, we have considered the f (R) Hu-Sawicki formulation for the dark energy component of the background Universe, in which the cosmic acceleration is driven by the no-Einsteinian geometrical terms if compared it to GR.
We have studied the dynamics in both cosmological scenarios to highlight the peculiarities of these models.To describe spherically symmetric deviations from homogeneity in the late Universe, we have used the separation of variables method to address the partial differential equations system for the first-order perturbation.Then, the radial component of such a reduction procedure has been analytically integrated, while the-time dependent part has required a numerical treatment.
The key difference between the two cosmological models studied in this work concerns the different form of the 0-1 component of the gravitational field equations.Indeed, in the ΛLTB model, this equation (19a) can be easily solved by removing one of the two free metric functions of the problem, hence the LTB metric takes the well-known simplified form (22). On the other hand, within the framework of the f (R) gravity in the Jordan frame, the 0-1 component given by Eq. (27a) intrinsically links the metric tensor components to the non-minimally coupled scalar field.Thus, the non-minimally coupling prevents the simplification above unlike GR, and we have to deal with two distinct metric functions α (t, r) and β (t, r) in the LTB metric given by Eq. (18).
Concerning the radial solutions, GR provides only a natural decay of the perturbations for large r values, following a power law, as it emerges from Eq. ( 53).Differently, in the f (R) modified gravity formulation, we obtain the Yukawa-like decaying for the radial perturbations of the scalar functions given by Eqs. ( 67), ( 68), (69), and (70).Furthermore, we have employed the f (R) HS model to describe the background Universe for comparison with the ΛCDM model, but we could have considered other viable f (R) modified gravity models.Indeed, we stress that our main result regarding the Yukawa-like radial perturbations does not depend on the f (R) functional form.
It should be noted that Yukawa-like radial profiles are recurring in f (R) theories, as evidenced, for instance, in other studies [115][116][117][118] by the presence of Yukawa-like corrections in the Newton potential with consequent implications for the dark matter problem and the flat rotation curves of galaxies.
In this paper, the different morphology of the radial solutions between the ΛLTB model and inhomogeneous f (R) cosmology must be regarded as the most relevant signature we fixed about the possibility to adopt the f (R) modified gravity scenario to describe the accelerating late Universe via spherically symmetric perturbations over a homogeneous background.In other words, the different radial profiles of local inhomogeneities in the Universe may be possible hints of a theory beyond GR.
For what concerned the time evolution of inhomogeneous perturbations both in the ΛLTB model and in the LTB solution as emerging from the f (R) gravity in the Jordan frame, the numerical analysis has allowed us to outline that it is always possible to obtain a non-divergent amplitude of the perturbations as time goes by, according to the reliable idea of a stable homogeneous and isotropic Universe in the near future.
Nevertheless, we are aware that our work has one limitation: the obtained radial solutions clearly diverge in the center of the LTB symmetry, where the observer, i.e., human location, is intended to be set.This feature simply suggests that our solution has a non-perturbative extension from a given large enough radial coordinate up to r = 0, which is an important task for future investigations in the late Universe.The present analysis fosters further studies about inhomogeneous cosmology when regarded as a local (non-linear) deformation of the FLRW geometry, as it may be also possibly pointed out by local measurements of the Hubble constant.
The present study may have a relevant impact on the observations of the large-scale structure of the Universe, when forthcoming missions, like the Euclid satellite [119], detect the clumpy galaxy distribution with greater accuracy.The comparison between the ΛCDM model and f (R) modified gravity theories through the investigation of spatial inhomogeneities may become a powerful tool to test the robustness of the cosmological concordance model.(B11c)

Declarations
We have tried to separate time-dependent and radial terms.However, this equations system exhibits a complicated structure, since it should be noted the occurrence of several mixed terms depending both on t and r.For instance, focusing on Eq. (B11a), terms on the left-hand side depend only on r, a mixed term is the first contribution on the right side, while the second contribution is only time-dependent.
Using again Eq. (B10), the linearized scalar field equation (61)  We noticed again the presence of mixed terms depending on t and r, which do not allow us to solve the equations using the separation of variables unless we rely on some simplifying assumptions.For instance, focusing on Eqs.(B11a) and (B13), if we require that the perturbations A p and B p follow a similar time evolution, i.e.
where λ 1 is a constant, then we are able to solve these two equations through the separation of variables method.Actually, we also write equivalently which is just Eq. (64a), since we can adjust constant term in the first-order perturbation theory.
Then, if we impose the assumption (B15) in Eq. (B11a), we obtain two equations, one in the variable t and the other one in r: which are just Eqs.(66a) and (65a), respectively.
In the same way, we can split Eq. (B13) into Eqs.(66e) and (65d) Ṗp + 3 ȧ ā P p + µ 2 ρ Ḃp = 0 , (B17a) (B18) At this point, noting a mixed term in the last contribution of the left-hand side, to proceed analytically with the separation of variables, we require an additional simplifying assumption, that is the proportionality between the radial evolution of the matter and scalar field perturbations: i.e., Eq. (64b), where λ 2 is the proportionality constant.
Hence, we can easily separate time and radial evolutions in Eq. (B18) to write the two differential equations (66d) and (65c), respectively:

Φp
which is just Eq. (66b).Finally, regarding the last equation of the system (B11) to be rewritten with the separation of variables, that is Eq.(B11c), if we combine it with Eqs.(B15), (B16a), (B16b), we obtain which are exactly Eqs.(66c) and (65b), respectively, where µ 4 is a constant.In conclusion, in this appendix, we have shown how the two simplifying assumptions (B15) and (B19) are suggested from the analysis of the equation system given by Eqs. ( 60), (61), and (62) to use the separation of variables.Finally, we have obtained separately two sets of equations ( 65) and (66) for the radial and time evolution of linear perturbations, which have been reported in Sect.6.2.Once we have split all field equations into time and space components, we can solve them to obtain linear order perturbations separately in Sect.6 have used the well-known relation ρ ∼ ā−3 for the matter component, the chain rule d dt = 1 t0 d dτ ≈ H 0 d dτ , and the usual definitions of the cosmological density parameters Ω m0 and Ω Λ0 .

Fig. 5
Fig. 5 Numerical solutions of linear perturbations in the Jordan frame of f (R) gravity in terms of the parameter τ .Top panel: the difference between the density contrast Ξ evaluated at a generic τ and τ = 1 today represented in units of 10 −11 and logarithmic scale.Bottom panel: the ratio between the linearly perturbed scalar field Φp (τ ) and the respective background scalar field φ (τ ) in units of 10 −13 .