1 Introduction

In recent decades, researchers came up with certain issues which cannot be expounded appropriately using general relativity (GR). A huge amount of unseen energy and matter famously known as Dark Energy (DE) and Dark Matter (DM) is found to be the main reason behind the issues arising in GR. The modified gravity (MG) theories endeavor to find possible description for the issues occurring in the standard cosmological model known as \(\Lambda \)CDM model by making alterations to GR. The validation of MG theories include the comparison of predictions of such theories with the already carried out solar system observations, e.g., accelerated cosmological expansion, rotation curves of galaxies and the cosmic microwave background spectrum. Since, we know that our current cosmological model, i.e., \(\Lambda \)CDM model embodies the Einstein’s field equations (EFEs) from GR as \(G_{\mu \nu }=\kappa T_{\mu \nu }\), where \(G_{\mu \nu }\) and \(T_{\mu \nu }\) symbolize the Einstein’s tensor and the energy-momentum tensor, respectively. The Einstein’s tensor incorporates the Riemann tensor, metric tensor and the Ricci curvature. The inclusion of these tensors provides a pathway to describe the distortion in spacetime arising due to matter and energy. Taking into account the current accelerated expansion, EFEs from GR seek insufficiency to elucidate the ongoing physical processes in our universe. The MG theories make use of an altered form of EFEs to predict the spacetime curvature including the effects of invisible DE and DM. Different MG theories exist with each one using a different criteria in order to provide a possible description of DM and DE. The f(R) gravity theory is one of the most famous MG theory which uses a straightforward approach to alter EFEs, i.e., the replacement of Ricci Scalar R by an arbitrary function f(R) of Ricci Scalar in the generalized Lagrangian of Einstein–Hilbert action.

Within most of the galaxies, there exist scattered clouds of dust. Due to the force of attraction between these dust particles, such clouds undergo a collapsing process and end up giving birth to a star. Over the course of time, this star changes its states depending upon the amount of mass it contains initially. This process is termed as stellar evolution. Scientists and researchers have been introducing novel ideas including building up a model to demonstrate the internal structure of a star. This has helped us in making predictions about the color, luminosity and the future evolutionary stage of a particular star. In order to construct such models, various techniques have been used, out of which structure scalars play a key role. These scalar functions emerges, when we split the Riemann tensor orthogonally and are known to exhibit rudimentary fluid properties like dissipative flux, pressure anisotropy and energy density inhomogeneity etc. In this way, such scalars help researchers in comprehending the formation and time-evolution of stars.

The deformity at large scale curvatures assists in the comprehension of current accelerated expansion of universe without utilizing DE. The f(R) theory is one of such deformed theories that opens a way to explain the dynamics at large scales. Sotiriou [1] inferred a relation between f(R) gravity and scalar-tensor theory using the Palatini approach. They inspected the circumstances under which both the theories become equivalent and investigated the consequences of such equivalence. Song et al. [2] discussed the evolution of linear perturbations in f(R) gravity models for accelerated expansion keeping in view the frame in which dynamical equations are of fourth order with minimally coupled data. Starobinsky [3] proposed a class of models producing a different and viable cosmology. They found that in the absence of matter, such models have de-Sitter and Minkowski spacetimes as particular solutions. A new issue regarding DE models based on f(R) gravity is also pointed out having a connection with scalarons (newly discovered massive particles). Capozziello et al. [4] explored the exact solutions for constant Ricci scalar and Ricci curvature scalar producing famous post-Minkowskian and post-Newtonian limits. They compared the results with the results of GR and also found first order solutions utilizing the perturbation scheme. Azadi et al. [5] discussed cylindrically symmetric vacuum solutions using Weyl coordinates keeping in view the metric f(R) gravity. They reduced modified Einstein’s equations into single equation exhibiting the fact that exact solutions can be constructed in correspondence with different f(R) gravity models. Capozziello et al. [6] showed the possibility of relating the cosmological parameters to the most recent values of f(R) along with its derivatives. They concluded that there is a possibility to relate the outcomes of the cosmography to the theoretical assumptions of f(R) cosmology. Bhatti and his collaborators [7,8,9,10,11,12] studied the effects of various modified gravity models on the dynamics as well as stability of different self-gravitating compact objects.

The anisotropic nature of pressure is all over the place in cosmos and much research has been carried out to comprehend its effects on the physical conduct of collapsing astrophysical objects. Raposo et al. [13] introduced the concept of studying the dynamical properties of self-gravitating anisotropic fluids by taking into account the covariant framework. They concluded that their work is specifically useful in blackhole paradigm. Bolokhov et al. [14] considered stationary cylindrically symmetric spacetimes along with anisotropic fluid distributions having directional pressures as the expected source of gravity. Keeping in view such sources, they investigated a way to obtain exact solutions by the splitting of Ricci tensor into rotational and static parts. Bowers and Liang [15] explored the significance of anisotropic equations of state by using the generalized form of hydrostatic equilibrium equation. They solved a numerical model having special form of anisotropy and concluded that few models lead to an increased value of surface redshift of same magnitude order as that of the fractional anisotropy.

Structure scalars are linked with the understanding of physical behavior of self-gravitating fluid distributions. To describe the evolution of astrophysical objects, many physicists utilized this concept. Herrera et al. [16] inspected the part played by the structure scalars for dissipative spherical fluid distributions in the presence of electromagnetic field. They considered neutral dust along with cosmological constant and discussed the physical outcomes of scalar functions. Herrera et al. [17] considered static axially symmetric sources and evaluated Einstein equations as well as the conservation equations having axial symmetry. They calculated the structure scalars for their case and determined the inhomogeneity factors and found few exact analytical solutions. Herrera et al. [18] calculated a set of stellar equations administering the composition and evolution of cylindrically symmetric self-gravitating dissipative fluids having anisotropic stresses. Such scalar quantities provide basic information, i.e., local anisotropy of pressure, energy density, dissipative flux, energy density inhomogeneity and active gravitational mass about the fluid under consideration. Recently, we have explore the effects of electromagnetic field on the structure formation of massive compact objects due to these structure scalars as well as on the polytropic models [19, 20].

The spherically symmetric solutions have drawn the attention of many scientists and physicists. Chodos and Detweiler [21] found the vacuum solution to five-dimensional Einstein’s equations and explored its significance keeping in view an already introduced Kaluza–Klein model. Wyman [22] reviewed that few scientists found spherically symmetric static solutions to a special kind of Einstein’s equations but they failed to evaluate explicit expressions for gravitational potentials. Multmäki and Vilja [23] explored the static spherically symmetric empty space solutions by taking into account the f(R) theory of gravity. They exhibited that for large class models, the Schwarzschild–de Sitter metric comes out to be an exact solution of Einstein’s stellar equations.

The phenomena in which an astrophysical object experiences certain changes with time is called stellar evolution. A star’s mass decides its fate, ranging from low to high, after facing evolutionary stages. Chabrier and Baraffe [24] investigated the formation and evolution of low mass stars and examined thermal and mechanical properties of such stars. They demonstrated the outcomes in the form of tables containing mass-radius-luminosity-effective temperature relations. Thorne [25] reformulated the relativistic equations of structure and evolution using a notation which is easily linked with the Newtonian theory. Lebreton [26] reviewed the outcomes of recent observations keeping in view the hippacros mission. They also presented expected and required theoretical advancements.

The part played by the physical variables on the time-evolution of variety of spacetimes has been discussed by many physicists. Kramer [27] explored a collapsing sphere model having heat flow in radial direction by taking the exterior solution to be Vaidya’s radiation field. Barreto [28] considered a radiating sphere and inspected the effects of viscosity on the gravitational collapse. It is concluded that such spherical model explodes viciously, when the viscosity is smaller and shear viscosity vanishes when the matter expands. Buchert [29] derived generalized Friedmann equations for spacetimes with irrotational dust and found a relationship between average scalar curvature and back reaction of such cosmology.

The active gravitational mass of any astrophysical object is termed as Tolman mass and it is the measure of the total energy budget of the astrophysical object. Devitt and Florides [30] derived modified Tolman mass-energy formula for time-independent and spherically symmetric systems. They found that when there is no surface of discontinuity, the evaluated formula exhibits significant properties. Abreu and Visser [31] demonstrated that a quasi-local Tolman mass in a volume can be reduced into Gauss-like surface integral for static spacetime. Herrera et al. [32] evaluated Tolman mass of a source after its equilibrium state and concluded that low departure from sphericity for compact objects result in increased or decreased value of Tolman mass depending on the equilibrium state of that system.

A set of congruences of worldlines must be specified in order to carry out observations. If the four-velocity of the fluid distribution is perpendicular to the group trajectories, then the system is said to be non-tilted, otherwise, it would be considered as tiled system. Much research has been carried out in order to explore the role of both the congruences. King and Ellis [33] investigated homogeneous cosmological models filled with perfect fluid and found relationship between well-known kinematical quantities (e.g., shear and rotation) and Bianchi classification of considered symmetry group. Triginer and Pavón [34] analyzed thermodynamical attributes of tilted fluids and evaluated the Gibbs equation as well as discussed the part played by the entropy. Barrow and Hervik [35] carried out the stability analysis of Bianchi type spatially homogeneous models incorporating tilted perfect fluids. Apostolopoulos [36] explored spatially homogeneous cosmologies of certain type incorporating tilted perfect fluid and demonstrated the outcomes by presenting few illustrative examples. Coley et al. [37] considered two time-like congruences and analyzed cosmological models having perfect fluid with a constant equation of state. They also discussed spatially and tilted spatially homogeneous models by taking into account the effects of cosmological constant.

This manuscript is devoted to comprehend the constitution and time-evolution of dissipative, self-gravitating spherically symmetric anisotropic stellar objects. Section 2 incorporates modified Einstein equations and kinematical quantities along with the expressions of Weyl tensor and the Riemann tensor. Two kinds of mass functions are evaluated and few stellar equations are also manipulated. In Sect. 3, the expressions for structure scalars in case of modified gravity theory are found by using the scheme of orthogonal splitting of Riemann tensor. Section 4 encompasses certain static stellar spheres as possible solutions to our line element. Section 5 concludes the findings of the manuscript.

2 Modified field equations

To arrive at the Einstein’s modified field equations, we alter the Einstein–Hilbert action as

$$\begin{aligned} S_{f(R)}=\frac{1}{2\kappa }\int d^4x \sqrt{-g}f(R)+S_M. \end{aligned}$$

It is evident that an arbitrary function f(R) of scalar curvature has taken place of the scalar curvature R utilized in standard action for GR. The alphabet g indicates the determinant of the metric tensor whereas \(\kappa \) symbolizes the coupling constant. The matter field action is denoted by \(S_M\). After some manipulations, we acquire modified Einstein’s gravitational equations as

$$\begin{aligned} G^\mu _\nu =\frac{\kappa }{f_R}( {T^\mu _\nu }^{(D)}+ T^\mu _\nu )\equiv {T^\mu _\nu }^{eff}. \end{aligned}$$
(1)

Bearing in mind the spherical configuration of astrophysical stellar object, we take into account dissipative and locally anisotropic fluid. The fluid is supposed to be contributing only in the interior region of the configuration. Using the Einstein frame of reference, we will work out the modified field equations. Moreover, it is assumed that the observer is not comoving with the fluid i.e. the reference frame of the observer is lorentz boosted in radial direction w.r.t the comoving one.

The line element for gravitating massive object is taken of the following form

$$\begin{aligned} ds^2= e^\nu dt^2-e^\lambda dr^2-r^2 d\theta ^2-r^2 sin^2 \theta d\phi ^2, \end{aligned}$$

where \(\nu \) and \(\lambda \) both are dependent on temporal coordinate t and radial coordinate r. The field equation yield the following independent components as

$$\begin{aligned} \kappa T_0^0&=f_R\left[ \frac{\lambda ' e^{-\lambda }}{r}+\frac{1}{r^2}-\frac{e^{-\lambda }}{r^2}\right] -\kappa \left. T^0_0\right. ^{(D)}, \end{aligned}$$
(2)
$$\begin{aligned} \kappa T_1^1&=f_R\left[ -\frac{\nu 'e^{-\lambda }}{r}+\frac{1}{r^2}-\frac{e^{-\lambda }}{r^2}\right] - \kappa \left. T^1_1\right. ^{(D)}, \end{aligned}$$
(3)
$$\begin{aligned} \kappa T_2^2&=f_R\left[ e^{-\nu }\left( \frac{\ddot{\lambda }}{2}+\frac{\dot{\lambda }^2}{4}-\frac{\dot{\lambda } \dot{\nu }}{4}\right) \right. \nonumber \\&\quad \left. +e^{-\lambda }\left( \frac{{\lambda }'}{2r}-\frac{{\nu }''}{2}-\frac{{\nu }'^2}{4}-\frac{{\nu }'}{2r}+\frac{{\lambda }' {\nu }'}{4}\right) \right] -\kappa \left. T^2_2\right. ^{(D)}, \end{aligned}$$
(4)
$$\begin{aligned} \kappa T_{10}&=\left[ \frac{\dot{\lambda }}{r}\right] f_R -\kappa \left. T_{10}\right. ^{(D)}. \end{aligned}$$
(5)

with dark source terms \(\left. T^0_0\right. ^{(D)},\left. T^1_1\right. ^{(D)},\left. T^2_2\right. ^{(D)}\) and \(\left. T_{10}\right. ^{(D)}\) being equal to

$$\begin{aligned} \left. T^0_0\right. ^{(D)}= & {} \frac{1}{\kappa }\left[ -\frac{\dot{\lambda } \dot{f_R} e^{-\nu }}{2}+f_R'' e^{-\lambda }\right. \\&\left. -\frac{\lambda ' f_R'e^{-\lambda }}{2}+\frac{2f_R'e^{-\lambda }}{r}+\frac{f-Rf_R}{2}\right] ,\\ \left. T^1_1\right. ^{(D)}= & {} \frac{1}{\kappa }\left[ \frac{\dot{\nu }\dot{f_R}e^{-\nu }}{2}+\ddot{f_R}e^{-\nu }+\frac{\nu ' f_R'e^{-\lambda }}{2}\right. \\&\left. +\frac{2e^{-\lambda }f_R'}{r}+\frac{(f-Rf_R)e^\lambda }{2}\right] ,\\ \left. T^2_2\right. ^{(D)}= & {} \frac{1}{\kappa }\left[ -\ddot{f_R}e^{-\nu }+\frac{f_R'e^{-\lambda }}{r}-\frac{\dot{\nu }\dot{f_R} e^{-\nu }}{2}\right. \\&\left. +\frac{\nu ' f_R'e^{-\lambda }}{2}+f_R'' e^{-\lambda }-\frac{\dot{\lambda } \dot{f_R} e^{-\nu }}{2}\right. \\&\left. -\frac{\lambda 'f_R'e^{-\lambda }}{2}+\frac{(f-Rf_R)}{2}\right] ,\\ \left. T_{10}\right. ^{(D)}= & {} \frac{1}{\kappa }\left[ \frac{\nu '\dot{f_R}}{2}+\frac{\dot{\lambda }f_R'}{2}-\dot{f_R'}\right] . \end{aligned}$$

Applying the Bondi formalism [38], the pure locally Minkowski coordinates symbolized as \((\tau ,x,y,z)\) are

$$\begin{aligned}&d\tau =e^{\nu /2}dt; \quad \quad dy=rd\theta ;\\&dx=e^{\lambda /2}dr; \quad \quad dz=r sin\theta d\phi . \end{aligned}$$

Utilizing the concept of tilted congruence, the Minkowskian components of the stress–energy tensor are expressed in the following form

$$\begin{aligned}&\bar{T}_0^0=T_0^0; \quad \quad \bar{T}_1^1=T_1^1;\\&\bar{T}_2^2=T_2^2; \quad \quad \bar{T}_3^3=T_3^3;\\&\bar{T}_{01}=e^{-(\nu +\lambda )/2}T_{01} . \end{aligned}$$

Taking into account a comoving observer, the components of stress–energy tensor in covariant form take the following form

$$\begin{aligned} \begin{pmatrix} \rho +\varepsilon &{}\quad -q-\varepsilon &{}\quad 0 &{}\quad 0\\ -q-\varepsilon &{}\quad P_r+\varepsilon &{}\quad 0 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad P_\bot &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 0 &{}\quad P_\bot \\ \end{pmatrix}. \end{aligned}$$

Application of Lorentz transformation to the above mentioned Minkowski components provide us with the following results

$$\begin{aligned} \bar{T}_0^0&=T_0^0=\frac{\rho +P_r \omega ^2}{1-\omega ^2}+\frac{2\omega q}{1-\omega ^2}+\frac{\varepsilon (1+\omega )}{1-\omega }, \end{aligned}$$
(6)
$$\begin{aligned} \bar{T}_1^1&=T_1^1=-\left( \frac{\rho \omega ^2+P_r}{1-\omega ^2}+\frac{2\omega q}{1-\omega ^2}+\frac{\varepsilon (1+\omega )}{1-\omega }\right) , \end{aligned}$$
(7)
$$\begin{aligned} \bar{T}_{01}&=e^{-(\nu +\lambda )/2}T_{01}\nonumber \\&=-\left( \frac{\omega e^{(\nu +\lambda )/2}(\rho +P_r)}{1-\omega ^2}+\frac{qe^{(\nu +\lambda )/2}(1+\omega ^2)}{1-\omega ^2}\right. \nonumber \\&\quad \left. +\frac{\varepsilon (1+\omega )e^{(\nu +\lambda )/2}}{1-\omega }\right) , \end{aligned}$$
(8)
$$\begin{aligned} \bar{T}_2^2&=T_2^2=-P_\bot ; \quad \bar{T}_3^3=T_3^3=-P_\bot . \end{aligned}$$
(9)

The time-derivative of radial coordinate i.e. dr/dt termed as coordinate velocity is interlinked with \(\omega \) by the following relation

$$\begin{aligned} \omega =\left( \frac{dr}{dt}\right) e^{\frac{(\lambda -\nu )}{2}}. \end{aligned}$$

Putting back Eqs. (6)–(9) into Eqs. (2)–(5), the following form of field equations is attained

$$\begin{aligned}&\frac{\rho +P_r \omega ^2}{1-\omega ^2}+\frac{2\omega q}{1-\omega ^2}+\frac{\varepsilon (1+\omega )}{1-\omega }\\&\quad =\frac{1}{\kappa }\left( \left[ \frac{\lambda ' e^{-\lambda }}{r}+\frac{1}{r^2}-\frac{e^{-\lambda }}{r^2}\right] f_R -\kappa \left. T^0_0\right. ^{(D)}\right) ,\\&\qquad -\left( \frac{\rho \omega ^2+P_r}{1-\omega ^2}+\frac{2\omega q}{1-\omega ^2}+\frac{\varepsilon (1+\omega )}{1-\omega }\right) \\&\quad =\frac{1}{\kappa }\left( \left[ -\frac{\nu ' e^{-\lambda }}{r}+\frac{1}{r^2}-\frac{e^{-\lambda }}{r^2}\right] f_R -\kappa \left. T^1_1\right. ^{(D)}\right) ,\\&\qquad -P_\bot =\frac{1}{\kappa }\left( \left[ e^{-\nu }\left( \frac{\ddot{\lambda }}{2}+\frac{\dot{\lambda }^2}{4}-\frac{\dot{\lambda } \dot{\nu }}{4}\right) \right. \right. \\&\qquad \left. \left. +e^{-\lambda }\left( \frac{{\lambda }'}{2r}-\frac{{\nu }''}{2}-\frac{{\nu }'^2}{4}-\frac{{\nu }'}{2r}+\frac{{\lambda }' {\nu }'}{4}\right) \right] f_R-\kappa \left. T^2_2\right. ^{(D)}\right) ,\\&\qquad -\left( \frac{\omega e^{(\nu +\lambda )/2}(\rho +P_r)}{1-\omega ^2}+\frac{qe^{(\nu +\lambda )/2}(1+\omega ^2)}{1-\omega ^2}\right. \\&\qquad \left. +\frac{\varepsilon (1+\omega )e^{(\nu +\lambda )/2}}{1-\omega }\right) = \frac{1}{\kappa }\left( \frac{\dot{\lambda }}{r}f_R-\kappa \left. T_{10}\right. ^{(D)}\right) . \end{aligned}$$

The four-velocity vector \(u^\beta \) symbolizing the velocity of the fluid inside the spherically symmetric configuration is taken to be as follows

$$\begin{aligned} u^\beta =\left( \frac{e^{-\nu /2}}{\sqrt{1-\omega ^2}},\frac{\omega e^{-\lambda /2}}{\sqrt{1-\omega ^2}},0,0\right) . \end{aligned}$$

Using four-velocity vector, we acquire four-acceleration components as follows

$$\begin{aligned} a_1\omega= & {} -a_0e^{(\lambda -\nu )/2}\\= & {} -\left( \frac{\omega ^2 \omega '}{(1-\omega ^2)^2}+\frac{\omega \nu '}{2(1-\omega ^2)^2}+ \frac{\omega ^2\dot{\lambda }e^{(\lambda -\nu )/2}}{2(1-\omega ^2)}\right. \\&\left. +\frac{\omega \dot{\omega }e^{(\lambda -\nu )/2}}{(1-\omega ^2)^2} \right) . \end{aligned}$$

For our case, shear tensor \(\sigma _{\beta \gamma }\) comes out to be as follows

$$\begin{aligned} \sigma _{\beta \gamma }=u_{\beta ;\gamma }+u_{\gamma ;\beta }-u_\beta a_\gamma -u_\gamma a_\beta -\frac{2\Theta h_{\beta \gamma }}{3}, \end{aligned}$$

with

$$\begin{aligned} h_{\beta \gamma }=g_{\beta \gamma }-u_\beta u_\gamma ; \quad \quad \Theta =u^\beta _{; \beta }. \end{aligned}$$

The expansion scalar is evaluated as follows

$$\begin{aligned} \Theta= & {} \frac{\dot{\lambda }e^{-\nu /2}}{2\sqrt{(1-\omega ^2)}}+ \frac{\omega \dot{\omega }e^{-\nu /2}}{(1-\omega ^2)^{3/2}}+ \frac{\omega \nu 'e^{-\lambda /2}}{2\sqrt{(1-\omega ^2)}}\\&+ \frac{\omega 'e^{-\lambda /2}}{(1-\omega ^2)^{3/2}}+ \frac{2\omega e^{-\lambda /2}}{r\sqrt{(1-\omega ^2)}}. \end{aligned}$$

We enlist the components of shear tensor as follows

$$\begin{aligned} \sigma _{00}&=\omega ^2 e^{\nu -\lambda } \sigma _{11},\\ \sigma _{01}&=-\omega e^{(\nu -\lambda )/2}\sigma _{11},\\ \sigma _{11}&=-\frac{2}{3(1-\omega ^2)^{3/2}}\left[ \dot{\lambda }e^\lambda e^{-\nu /2}+\frac{2\omega \dot{\omega } e^\lambda e^{-\nu /2}}{(1-\omega ^2)}\right. \\&\quad \left. +\omega \nu 'e^{\lambda /2}+\frac{2\omega 'e^{\lambda /2}}{(1-\omega ^2)}-\frac{2\omega e^{\lambda /2}}{r} \right] , \\ \sigma _{22}&=-\frac{r^2e^{-\lambda }(1-\omega ^2)}{2} \sigma _{11},\\ \sigma _{33}&=-\frac{r^2e^{-\lambda }(1-\omega ^2)}{2} sin^2\theta \sigma _{11}. \end{aligned}$$

Shear tensor can also be expressed by using an alternative form as

$$\begin{aligned} \sigma _{\beta \gamma }=\frac{\sigma }{2} \left( s_\beta s_\gamma +\frac{1}{3}h_{\beta \gamma } \right) , \end{aligned}$$

with \(\sigma \) being equal to the following value

$$\begin{aligned} \sigma= & {} -\frac{1}{\sqrt{(1-\omega ^2)}}\left[ \dot{\lambda }e^{-\nu /2}+ \frac{2\omega \dot{\omega }e^{-\nu /2}}{(1-\omega ^2)}\right. \\&\left. +\omega \nu ' e^{-\lambda /2}+\frac{2\omega ' e^{-\lambda /2}}{(1-\omega ^2)}-\frac{2\omega e^{-\lambda /2}}{r} \right] . \end{aligned}$$

We define another four-vector \(s^\beta \) as follows

$$\begin{aligned} s^\beta =\left( \frac{\omega e^{-\nu /2}}{\sqrt{(1-\omega ^2)}},\frac{e^{-\lambda /2}}{\sqrt{(1-\omega ^2)}},0,0\right) . \end{aligned}$$

It is found that \(s^\beta \) satisfies the following equations

$$\begin{aligned} s^\beta u_\beta =0; \quad \quad s^\beta s_\beta =-1. \end{aligned}$$

The stress–energy tensor which depicts the characteristics of the fluid distribution is stated as follows

$$\begin{aligned} T^\beta _\gamma =\tilde{\rho }u^\beta u_\gamma -\check{P}h^\beta _\gamma +\Pi ^\beta _\gamma + \tilde{q}(s^\beta u_\gamma +s_\gamma u^\beta ), \end{aligned}$$

with

$$\begin{aligned} \Pi ^\beta _\gamma&=\Pi \left( s^\beta s_\gamma +\frac{1}{3}h^\beta _\gamma \right) ; \quad \check{P}=\frac{\tilde{P}_r +2P_\bot }{3};\\ \tilde{q}^\beta&=\tilde{q}s^\beta ; \quad \tilde{\rho }=\rho +\varepsilon ; \quad \tilde{P_r}=P_r+\varepsilon ;\\ \tilde{q}&=q+\varepsilon ; \quad \Pi =\tilde{P}_r-P_\bot . \end{aligned}$$

To manipulate the matching conditions, Vaidya spacetime is taken to be the exterior metric and its expression is stated as

$$\begin{aligned} ds^2=\left( \frac{1-2M(v)}{\mathfrak {R}}\right) dv^2 +2dvd\mathfrak {R}-\mathfrak {R}^2(d\theta ^2+ sin^2 \theta d\phi ^2), \end{aligned}$$

here v and \(\mathfrak {R}\) are symbols for retarded time coordinate and null coordinate, respectively. To attain smooth matching of the interior and exterior spacetimes over the boundary surface \(r=r_\Sigma \), the continuity of first and second fundamental forms are kept in view and the following outcomes are produced.

$$\begin{aligned}&e^{\nu _\Sigma }=\left( 1-\frac{2M}{R}\right) _{r_\Sigma }; e^{\lambda _\Sigma }=\left( 1-\frac{2M}{R}\right) _{r_\Sigma },\nonumber \\&\quad (P_r^{(eff)})_\Sigma = (q^{(eff)})_\Sigma , \end{aligned}$$
(10)

with \(\Sigma \) signifying the fact that the quantity is calculated on the boundary surface. Here, we have \(P_r^{(eff)}=\kappa P_r+\tau ^1_1\) and \(q^{(eff)}=\kappa q+\tau _{10}\). The quantities \(\tau ^1_1\) and \(\tau _{10}\) defined below are not tensors rather they are just notations used to express the contribution of dark source terms in the respective stellar equation.

$$\begin{aligned} \tau ^1_1&=-\frac{\dot{\nu }\dot{f_R}e^{-\nu }}{2}-\ddot{f_R}e^{-\nu }-\frac{\nu ' f_R'e^{-\lambda }}{2}\\&\quad -\frac{2e^{-\lambda }f_R'}{r}-\frac{(f-Rf_R)e^\lambda }{2},\\ \tau _{10}&=-\dot{f}_R'+\frac{\nu '\dot{f}_R}{2}+\frac{\dot{\lambda }f_R'}{2}. \end{aligned}$$

2.1 The Weyl tensor and the Riemann tensor

Besides measuring the curvature of spacetime, the Weyl tensor \(C^\rho _{\alpha \beta \mu }\) can also be utilized in an alternative definition of Riemann tensor along with Ricci tensor \(R_{\alpha \beta }\) and scalar curvature R as

$$\begin{aligned} R^\rho _{\eta \pi \mu }= & {} C^\rho _{\eta \pi \mu }+\frac{1}{2}R^\rho _\pi g_{\eta \mu }+\frac{1}{2}R_{\eta \pi }\delta ^\rho _\mu +\frac{1}{2}R_{\eta \mu }\delta ^\rho _\pi \nonumber \\&-\,\frac{1}{2}R^\rho _\mu g_{\eta \pi }-\frac{1}{6}R\left( \delta ^\rho _\pi g_{\eta \mu }-g_{\eta \pi }\delta ^\rho _\mu \right) . \end{aligned}$$
(11)

Since, our line element is spherically symmetric, the magnetic part of the Weyl tensor comes out to be equal to zero. So, we utilize only its electric part given as follows for further discussion, i.e., \(E_{\sigma \beta }=C_{\sigma \beta \gamma \delta }u^\gamma u^\delta \) as

$$\begin{aligned} C_{\xi \nu \kappa \lambda }=E^{\beta \delta }u^\alpha u^\gamma (g_{\xi \nu \alpha \beta } g_{\kappa \lambda \gamma \delta }-\eta _{\xi \nu \alpha \beta } \eta _{\kappa \lambda \gamma \delta }), \end{aligned}$$

where \(g_{\xi \nu \alpha \beta }\) equals \(g_{\xi \nu \alpha \beta }=g_{\xi \alpha }g_{\nu \beta }-g_{\xi \beta }g_{\nu \alpha }\) and \(\eta _{\mu \nu \alpha \beta }\) signifies well-known Levi-Civita tensor. \(E_{\sigma \beta }\) has the following alternative form

$$\begin{aligned} E_{\sigma \beta }=E\left( s_\sigma s_\beta +\frac{1}{3} h_{\sigma \beta } \right) , \end{aligned}$$
(12)

with E being equal to

$$\begin{aligned} E&=\frac{\ddot{\lambda }e^{-\nu }}{4}+\frac{\dot{\lambda }^2 e^{-\nu }}{8}-\frac{\dot{\lambda }\dot{\nu } e^{-\nu }}{8}-\frac{{\nu }''e^{-\lambda }}{4}-\frac{\nu '^2 e^{-\lambda }}{8}+\frac{\lambda '\nu 'e^{-\lambda }}{8}\\&\quad +\,\frac{\nu 'e^{-\lambda }}{4r}- \frac{\lambda 'e^{-\lambda }}{4r}-\frac{e^{-\lambda }}{2r^2}+\frac{1}{2r^2}. \end{aligned}$$

Observations reveal that E satisfies the following equations

$$\begin{aligned} E^\beta _\beta =0; \quad \quad E_{\eta \pi }=E_{(\eta \pi )}; \quad \quad E_{\eta \pi }u^\pi =0. \end{aligned}$$

2.2 The mass function, Tolman mass and structural equations

In this section, we will explore the mass function formulated by Tolman for our massive object under the influence of modified gravity. Also, we will formulate the evolution equations for shear and expansion scalars in terms of structure scalars.

2.2.1 The mass function

The mass function m exhibiting the spherical geometry is known to be equal to

$$\begin{aligned} 1-e^{-\lambda }=\frac{2m}{r}. \end{aligned}$$
(13)

Making use of the modified field equations along with Eqs. (11) and (12) we deduce

$$\begin{aligned} \frac{3m}{r^3}=\frac{\kappa }{2f_R}\left( \tilde{\rho }-\tilde{P_r}+\tilde{P_\bot }\right) +E-\left( \frac{\tau ^0_0+\tau ^1_1}{2f_R}\right) . \end{aligned}$$
(14)

Here, we have

$$\begin{aligned} \tau ^0_0=\frac{\dot{\lambda } \dot{f_R} e^{-\nu }}{2}-f_R'' e^{-\lambda }+\frac{\lambda ' f_R'e^{-\lambda }}{2}-\frac{2f_R'e^{-\lambda }}{r}-\frac{f-Rf_R}{2}. \end{aligned}$$

Another definition of mass function m(r) is evaluated as

$$\begin{aligned} m(r)=\frac{\kappa }{2} \int ^r_0 \frac{r^2}{f_R} \left( T^0_0-\frac{\tau ^0_0}{\kappa }\right) dr. \end{aligned}$$
(15)

This equation links the static mass of interior spherical geometry with the energy density in the background of modified theory of gravity. Using the modified Einstein’s equations along with Eqs. (11), (12) and (15)

$$\begin{aligned} m= & {} \frac{\kappa r^3 }{6f_R}\left( T^0_0+T^1_1-T^2_2\right) \nonumber \\&+\,\frac{r^3E}{3}-\frac{r^3(\tau ^0_0+\tau ^1_1)}{6f_R}. \end{aligned}$$
(16)

Our chief target is to acquire an expression of tidal forces E in terms of the energy density inhomogeneity \({T^0_0}'\) in addition with pressure anisotropy \(T^2_2-T^1_1\) and dark source ingredients \(\tau ^0_0,\tau ^1_1\). For this purpose, we differentiate Eq. (16) and then integrate the obtained result. Thus, taking derivative w.r.t. the radial coordinate r, we acquire

$$\begin{aligned} \left( \frac{r^3 E}{3}\right) '= & {} -\frac{\kappa r^3}{6f_R}\left[ (T^0_0)'+\frac{3\tau ^0_0}{\kappa r}\right] \nonumber \\&+\left[ \frac{\kappa r^3}{6f_R}(T^2_2-T^1_1)\right] '+\left[ \frac{r^3(\tau ^0_0+\tau ^1_1)}{6f_R}\right] '. \end{aligned}$$
(17)

Now, integration of above equation provides us the following outcome

$$\begin{aligned} E= & {} -\frac{\kappa }{2r^3}\int _0^r \frac{r^3}{f_R}\left[ (T^0_0)'+\frac{3\tau ^0_0}{\kappa r}\right] dr\nonumber \\&+\,\frac{\kappa }{2f_R}(T^2_2-T^1_1)+\frac{(\tau ^0_0+\tau ^1_1)}{2f_R}. \end{aligned}$$
(18)

The mass function takes the following value

$$\begin{aligned} m(r,t)= & {} \frac{\kappa r^3 T_0^0}{f_R}-\frac{\kappa }{6} \int _0^r\frac{r^3}{f_R}\left[ (T^0_0)'+\frac{3\tau ^0_0}{\kappa r}\right] dr\nonumber \\&+\,\frac{r^3(\tau ^0_0+\tau ^1_1)}{6f_R}. \end{aligned}$$
(19)

Unlike Eq. (15), this equation signifies an important relationship between non-static interior mass m(r, t) of spherically symmetric configuration and rudimentary fluid properties \(T_0^0\), \({T_0^0}'\). Based on the following criteria, we can express stress–energy tensor components as \(T^0_0=\tilde{\rho }\) and \(T^1_1=\Pi \).

  1. (i)

    when \(\omega \) and time derivatives of all quantities vanish i.e. the static regime.

  2. (ii)

    If \(\omega ^2\approx \dot{\omega }\approx \ddot{\nu }\approx \ddot{\lambda }\approx \dot{\nu }\dot{\lambda } \approx \dot{\lambda }^2\approx 0\) i.e. the quasistatic regime.

  3. (iii)

    When \(\omega \approx \dot{\nu }\approx \dot{\lambda }\approx 0\) but \(\dot{\omega }\ne 0\)i.e. just after the departure of system from its state of equilibrium.

    Taking into account these three cases, Eqs. (18) can be re-written as follows

    $$\begin{aligned} E=-\frac{\kappa }{2r^3} \int ^r_0\frac{r^3}{f_R}\left[ \tilde{\rho }'+\frac{3 T_0^0}{f_R}\right] dr+\frac{\kappa \Pi }{2f_R}+\left( \frac{\tau ^0_0+\tau ^1_1}{2f_R}\right) , \end{aligned}$$
    (20)

    and Eq. (19) provides the following result

    $$\begin{aligned} m(r,t)=\frac{\kappa r^3 \tilde{\rho }}{f_R}-\frac{\kappa }{6} \int _0^r\frac{r^3}{f_R}\left[ (\tilde{\rho })'+\frac{3\tau ^0_0}{\kappa r}\right] dr+\frac{r^3(\tau ^0_0+\tau ^1_1)}{6f_R}. \end{aligned}$$
    (21)

    It can be seen clearly that Eq.(20) links Weyl tensor with the physical characteristics of fluid, e.g., energy density inhomogeneity and pressure anisotropy along with extra curvature terms. Equation (21) exhibits that mass function can be written in terms of homogeneous energy density in the background of f(R) theory of gravity.

2.2.2 The Tolman mass

To learn about the total energy budget of a celestial configuration, Tolman mass comes in handy. For spherical configurations, its expression is stated as

$$\begin{aligned} m_T&=\frac{\kappa }{2}\int ^{r_\Sigma }_0 r^2 e^{(\nu +\lambda )/2}\left( \left. T^0_0\right. ^{eff}-\left. T^1_1\right. ^{eff}-2\left. T^2_2\right. ^{eff}\right) dr\\&\quad +\frac{1}{2}\int ^{r_\Sigma }_0 f_R(R_0) r^2 e^{(\nu +\lambda )/2}\\&\quad \times \frac{\partial }{\partial t}\left[ \frac{\partial L}{\frac{\partial }{\partial t}[\partial (g^{\alpha \beta } \sqrt{-g})]}\right] g^{\alpha \beta }dr. \end{aligned}$$

with L and \(R_0\) symbolizing the gravitational Lagrangian density and the constant scalar curvature respectively. The expression for the mass function present inside sphere of radius r becomes

$$\begin{aligned} m_T&=\frac{\kappa }{2}\int ^r_0 r^2 e^{(\nu +\lambda )/2}(\left. T^0_0\right. ^{eff}-\left. T^1_1\right. ^{eff}-2\left. T^2_2\right. ^{eff})dr\\&\quad +\,\frac{1}{2}\int ^r_0 f_R(R_0) r^2 e^{(\nu +\lambda )/2}\\&\quad \times \, \frac{\partial }{\partial t}\left[ \frac{\partial L}{\frac{\partial }{\partial t}[\partial (g^{\alpha \beta } \sqrt{-g})]}\right] g^{\alpha \beta }dr, \end{aligned}$$

which in turn produces the following outcome

$$\begin{aligned} m_T= & {} e^{(\nu +\lambda )/2}\left[ m(r,t)f_R-\frac{\kappa r^3}{2}\left( T^1_1+\left. T^1_1\right. ^{(D)}\right) \right] \nonumber \\&-\,\frac{1}{2} \int _0^r r^2e^{(\lambda -\nu )/2} \dot{\lambda } \dot{f_R} dr. \end{aligned}$$
(22)

Substituting the values for \(T^1_1\) and m, we attain

$$\begin{aligned} m_T=\frac{\nu 'r^2f_R e^{(\nu -\lambda )/2}}{2}-\frac{1}{2} \int _0^r r^2 e^{(\lambda -\nu )/2}\dot{\lambda } \dot{f_R} dr. \end{aligned}$$

In static gravitational field, the instantaneous acceleration of a test particle is given by

$$\begin{aligned} a=\frac{m_Te^{-\nu /2}}{r^2}. \end{aligned}$$

We acquire another suitable expression for Tolman mass \(m_T\) as follows

$$\begin{aligned} rm_T'-3m_T&=\frac{r^3e^{(\nu -\lambda )/2}}{2}\left( \nu 'f_R'-\dot{\lambda } \dot{f_R}\right) +r^3f_R e^{(\lambda -\nu )/2}\nonumber \\&\quad \times \left[ \frac{\ddot{\lambda }}{2}+\frac{\dot{\lambda }^2}{4}-\frac{\dot{\lambda }\dot{\nu }}{4}\right] +\frac{r^3e^{(\lambda +\nu )/2}}{2}\nonumber \\&\quad \times \left[ \kappa \left( T^1_1+\left. T^1_1\right. ^{(D)}-T^2_2-2\left. T^2_2\right. ^{(D)}\right) -Ef_R\right] \nonumber \\&\quad +\frac{3}{2} \int _0^r r^2\dot{\lambda } \dot{f_R} e^{(\lambda -\nu )/2} dr. \end{aligned}$$
(23)

Applying integration, we deduce

Substituting the value for E in the above equation, we get

The components of stress–energy tensor are \(T^0_0=\tilde{\rho }\) and \(T^2_2-T^1_1=\Pi \) by keeping in view the three regimes stated before Eq. (20). The last equation exhibits that Tolman mass can be linked with its value for the case of isotropic fluid and homogeneous energy density.

2.2.3 Structure and evolution Stellar equations

This section enlists certain stellar equations that assist in the comprehension of formation, structure and time-evolution of celestial objects. In the coming sections, these equations are further written in terms of scalar functions to gain insight into the fundamental characteristics of fluid.

It is worth noticing that the operators \(*\) and \(\dagger \) are just the notations used to denote the derivatives of the quantities with a slight difference as follows

$$\begin{aligned} f^\dag = f_{,\mu }u^\mu ,\quad \quad \quad f^* = f_{,\mu }s^\mu . \end{aligned}$$

Firstly, we evaluate the first conservation equation in the context of f(R) theory as

$$\begin{aligned}&\tilde{\rho }^*+(\tilde{\rho }+\tilde{P}_r)\theta -\frac{2}{3}\left( \theta +\frac{\sigma }{2}\right) \Pi + \tilde{q}^\dagger +2\tilde{q}\nonumber \\&\quad \times \left( a+\frac{s^1}{r}\right) +\chi _1=0, \end{aligned}$$
(24)

with \(\chi _1\) being equal to

$$\begin{aligned} \chi _1&=\dot{\left. T^0_0\right. ^{(D)}}+\left( \frac{(\nu '-\lambda ')e^{(\nu -\lambda )/2}}{2(1-\omega ^2)}+\frac{2\omega \omega 'e^{(\nu -\lambda )/2}}{(1-\omega ^2)^2}\right) \\&\quad \left. T^1_0\right. ^{(D)}+ \acute{\left. T^1_0\right. ^{(D)}}+\frac{\dot{\lambda }\left. T^0_0\right. ^{(D)}}{2}+\frac{{\lambda }'\left. T^1_0\right. ^{(D)}}{2}\\&\quad +\frac{2\left. T^1_0\right. ^{(D)}}{r} +\frac{\nu 'e^{\nu -\lambda }\left. T^0_1\right. ^{(D)}}{2}-\frac{\dot{\lambda }\left. T^1_1\right. ^{(D)}}{2}+ \frac{\dot{f_R}}{f_R}\\&\quad \times \left[ \frac{\tilde{\rho }+\check{P}\omega ^2}{1-\omega ^2}+\frac{2\Pi \omega ^2}{3}+\frac{2\omega \tilde{q}}{1-\omega ^2}+\left. T^0_0\right. ^{(D)}\right] + \frac{f_R'}{f_R}\\&\quad \times \left[ \frac{e^{(\nu -\lambda )/2}}{1-\omega ^2}\left( (\tilde{\rho }+\check{P})\omega +\frac{2\Pi \omega }{3} +\tilde{q}(1+\omega ^2)\right) +\left. T^1_0\right. ^{(D)}\right] , \end{aligned}$$

while

$$\begin{aligned} \left. T^1_0\right. ^{(D)}=\frac{1}{\kappa }\left( -\frac{\nu '\dot{f_R}e^{-\nu }}{2}-\frac{\dot{\lambda }f_R'e^{-\nu }}{2}+e^{-\nu }\dot{f_R'}.\right) . \end{aligned}$$

The second conservation equation comes out to be as follows

$$\begin{aligned}&\tilde{P}_r^\dagger +(\tilde{\rho }+\tilde{P}_r)a+\frac{2s^1\Pi }{r}\nonumber \\&\quad -\frac{\tilde{q}}{3}(\sigma -4\theta )+\tilde{q}^*+\chi _2=0, \end{aligned}$$
(25)

where \(\chi _2\) has the following value

$$\begin{aligned} \chi _2&=\left. \dot{T^0_1}\right. ^{(D)}+\left. \acute{T^1_1}\right. ^{(D)}+\left( \frac{\nu '}{2}+\frac{2}{r}\right) \left. T^1_1\right. ^{(D)}\\&\quad -\frac{\nu '\left. T^0_0\right. ^{(D)}}{2}+\frac{\dot{\lambda }e^{\lambda -\nu }\left. T^1_0\right. ^{(D)}}{2(1-\omega ^2)}-\frac{2\left. T^2_2\right. ^{(D)}}{r}\\&\quad +\frac{\dot{\nu }\left. T^0_1\right. ^{(D)}}{2}+\frac{\dot{f_R}}{f_R}\left[ \frac{e^{(\lambda -\nu )/2}}{1-\omega ^2}\left( (\tilde{\rho }+\check{P})\omega \right. \right. \\&\quad \left. \left. +\frac{2\Pi \omega }{3}+\tilde{q}(1+\omega ^2)\right) +\left. T^0_1\right. ^{(D)}\right] \\&\quad +\frac{f_R'}{f_R}\left[ \frac{\tilde{\rho }\omega ^2}{1-\omega ^2} +\frac{\check{P}}{1-\omega ^2}+\frac{2\Pi }{3(1-\omega ^2)}\right. \\&\quad \left. +\frac{2\omega \tilde{q}}{1-\omega ^2}+\left. T^1_1\right. ^{(D)}\right] . \end{aligned}$$

The Raychaudhuri equation used for describing the motion of matter particles is stated in our case as

$$\begin{aligned}&\theta ^*+\frac{\theta ^2}{3}+\frac{\sigma ^2}{6}-a^\dagger -a^2-\frac{2as^1}{r}\nonumber \\&\quad =-\frac{\kappa }{2}(\tilde{\rho }+3\tilde{P}_r)+k\Pi +\chi _3, \end{aligned}$$
(26)

with the value of \(\chi _3\) as

$$\begin{aligned} \chi _3= & {} \left[ \frac{f-Rf_R}{2}+\frac{\Box f_R}{2}-\frac{u^\alpha u^\delta \nabla _\alpha \nabla _\delta f_R}{2}\right. \\&\left. -\frac{u^\beta u_\theta \nabla ^\theta \nabla _\beta f_R}{2}+2u_\theta u^\delta \nabla ^\theta \nabla _\delta f_R\right] . \end{aligned}$$

Now, using Ricci identities, we deduce the following result

$$\begin{aligned} \left( \sigma +\frac{\theta }{2}\right) ^\dagger =-\frac{3\sigma s^1}{2r}+\frac{3}{2}\sqrt{\left( \kappa ^2 \tilde{q}^2+2\psi _1\right) }, \end{aligned}$$
(27)

where \(\psi _1\) takes the following value

$$\begin{aligned} \psi _1=\kappa \tilde{q}u^\delta s_\rho \nabla ^\rho \nabla _\delta f_R-\frac{(\nabla ^\rho \nabla _\delta f_R)^2}{2}h_{\rho \rho } u^\delta u^\delta . \end{aligned}$$

Further, we make use of modified field equations along with Eq. (11) and Ricci identities to produce the following outcome

$$\begin{aligned}&-a^\dagger -a^2-\frac{\sigma ^*}{2}-\frac{\theta \sigma }{2}+\frac{as^1}{r}+\frac{\sigma ^2}{12}= E+\frac{\kappa \Pi }{2}\nonumber \\&\quad +\frac{\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }, \end{aligned}$$
(28)

where

$$\begin{aligned} \zeta _{\sigma \pi }=\frac{h^\mu _\sigma h^\nu _\pi \nabla _\mu \nabla _\nu f_R}{2}-\frac{h_{\sigma \pi }u_\theta u^\delta \nabla ^\theta \nabla _\delta f_R}{6}-\frac{h_{\sigma \pi }\Box f_R}{6}. \end{aligned}$$

Lastly, we use Bianchi identities in terms of the Weyl tensor to get the following two results

$$\begin{aligned}&\left( \frac{\kappa P_r}{2f_R}+\frac{3m}{r^3}\right) \left( \theta +\frac{\sigma }{2}\right) \nonumber \\&\qquad +\left[ \left( E-\frac{\kappa \Pi }{2}-\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right) +\left( \frac{\kappa \tilde{\rho }}{2}+\frac{\chi _4}{2}\right) \right] ^*\nonumber \\&\quad =-\frac{3s^1}{2r}\left( \sqrt{\left( \kappa ^2 \tilde{q}^2+2\psi _1\right) }\right) , \end{aligned}$$
(29)

with

$$\begin{aligned} \xi _{\sigma \pi }&=\frac{\epsilon ^{\epsilon \rho }_\sigma \epsilon _{\pi \mu \rho \nabla ^\mu \nabla _\epsilon f_R}}{8}+\frac{\epsilon ^{\epsilon \rho }_\sigma \epsilon _{\mu \epsilon \pi \nabla ^\mu \nabla _\rho f_R}}{8}\\&\quad +\frac{\epsilon ^{\epsilon \rho }_\sigma \epsilon _{\theta \pi \rho \nabla ^\theta \nabla _\epsilon f_R}}{8}+\frac{h_{\sigma \pi }h^\epsilon _\mu \nabla ^\mu \nabla _\epsilon f_R}{12}\\ {}&\quad +\frac{h_{\sigma \pi }h^\epsilon _\theta \nabla ^\theta \nabla _\epsilon f_R}{12}+\frac{h_{\sigma \pi }h^\rho _\mu \nabla ^\mu \nabla _\rho f_R}{12}\\ {}&\quad +\frac{h_{\sigma \pi }h^\rho _\theta \nabla ^\theta \nabla _\rho f_R}{12}, \end{aligned}$$

and

$$\begin{aligned}&\left[ \left( \frac{\kappa \tilde{\rho }}{2}+\frac{\chi _4}{2}\right) +\left( E-\frac{\kappa \Pi }{2}-\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right) \right] ^\dagger \nonumber \\&\quad =\frac{3s^1}{r}\left( \frac{\kappa \Pi }{2}-E+\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right) \nonumber \\&\qquad +\left( \theta +\frac{\sigma }{2}\right) \sqrt{\left( \frac{\kappa ^2\tilde{q}^2}{2} +\psi _1\right) }, \end{aligned}$$
(30)

where we have

$$\begin{aligned} \chi _4= & {} \left[ -\frac{f-Rf_R}{2}-\frac{h_\pi ^\epsilon \nabla ^\pi \nabla _\epsilon f_R}{4}-\frac{h^\rho _\pi \nabla ^\pi \nabla _\rho f_R}{4}\right. \\&\left. -\frac{h^\epsilon _\theta \nabla ^\theta \nabla _\epsilon f_R}{4}-\frac{h^\rho _\theta \nabla ^\theta \nabla _\rho f_R}{4}\right] . \end{aligned}$$

Equation (14) can also be written as

$$\begin{aligned} \frac{3m}{r^3}=\frac{1}{f_R}\left[ \frac{\kappa \tilde{\rho }}{2}+\frac{\kappa }{2}(P_\bot -\tilde{P}_r) \right] +E-\left( \frac{\tau ^0_0+\tau ^1_1}{2f_R}\right) , \end{aligned}$$
(31)

where we have \(a^\mu = as^\mu \).

3 The orthogonal breakdown of the Riemann tensor

This section deals with the evaluation of certain scalar functions by using the formalism in which the Riemann tensor is split orthogonally. For this purpose, we take into account the following tensor quantities

$$\begin{aligned}&Y_{\sigma \pi }=R_{\sigma \gamma \pi \delta }u^\gamma u^\delta ,\\&Z_{\sigma \pi }=^*R_{\sigma \gamma \pi \delta }u^\gamma u^\delta =\frac{1}{2}\eta _{\sigma \gamma \epsilon \rho }R^{\epsilon \rho }_{\pi \delta }u^\gamma u^\delta ,\\&X_{\sigma \pi }=^*R^*_{\sigma \gamma \pi \delta }u^\gamma u^\delta =\frac{1}{2}\eta _{\sigma \gamma }^{\epsilon \rho }R^*_{\epsilon \rho \pi \delta }u^\gamma u^\delta , \end{aligned}$$

with

$$\begin{aligned} R^*_{\sigma \pi \gamma \delta }=\frac{1}{2}\eta _{\epsilon \rho \gamma \delta }R^{\epsilon \rho }_{\sigma \pi }. \end{aligned}$$

Using the modified Einstein’s equations, we acquire

$$\begin{aligned} R^{\sigma \gamma }_{\pi \delta }&=C^{\sigma \gamma }_{\pi \delta }+2\kappa \left. T^{(eff)}\right. ^{[\sigma }_{[\pi }\delta ^{\gamma ]}_{\delta ]}\nonumber \\&\quad +\kappa \left. T^{(eff)}\right. \left( \frac{1}{3}\delta ^\sigma _{[\pi }\delta ^\gamma _{\delta ]}-\delta ^{[\sigma }_{[\pi }\delta ^{\gamma ]}_{\delta ]} \right) . \end{aligned}$$
(32)

Utilizing Eq. (32) and stress–energy tensor components, we perform the splitting of Riemann tensor into three parts as follows

$$\begin{aligned} R^{\sigma \gamma }_{\pi \delta }=R^{\sigma \gamma }_{(I)\pi \delta }+R^{\sigma \gamma }_{(II)\pi \delta }+R^{\sigma \gamma }_{(III)\pi \delta }, \end{aligned}$$

with each part defined as

All the three parts satisfy

$$\begin{aligned} \epsilon _{\sigma \gamma \pi }=u^\mu \eta _{\mu \sigma \gamma \pi }; \quad \quad \quad \epsilon _{\sigma \gamma \pi }u^\pi =0. \end{aligned}$$

Other important relations include

$$\begin{aligned}&\epsilon ^{\mu \gamma \nu }\epsilon _{\nu \sigma \pi }=u^\theta u^\rho \eta ^{\mu \gamma \nu }_{\rho }\eta _{\theta \nu \sigma \pi },\\&\epsilon ^{\mu \gamma \nu }\epsilon _{\nu \sigma \pi }=\delta ^\gamma _\sigma h^\mu _\pi -\delta ^\mu _\sigma h^\gamma _\pi +u_\sigma \left( u^\mu \delta ^\gamma _\pi -\delta ^\mu _\pi u^\gamma \right) . \end{aligned}$$

Upon contraction of \(\mu \) with \(\sigma \), we attain

$$\begin{aligned} \epsilon ^{\mu \gamma \nu }\epsilon _{\nu \mu \beta }=-2h^\gamma _\beta . \end{aligned}$$

We procure the following explicit expressions for the three tensor quantities \(X_{\sigma \pi },~Y_{\sigma \pi }\) and \(Z_{\sigma \pi }\)

$$\begin{aligned} X_{\sigma \pi }&=\frac{\kappa \tilde{\rho }}{3}h_{\sigma \pi }+\frac{\kappa \Pi _{\sigma \pi }}{2}-E_{\sigma \pi }+\psi _{\sigma \pi }, \end{aligned}$$
(33)
$$\begin{aligned} Y_{\sigma \pi }&=\frac{\kappa \tilde{\rho }}{6}h_{\sigma \pi }+\frac{\kappa \check{P}}{2}h_{\sigma \pi } +E_{\sigma \pi }+\frac{\kappa \Pi _{\sigma \pi }}{2}+\varphi _{\sigma \pi }, \end{aligned}$$
(34)
$$\begin{aligned} Z_{\sigma \pi }&=\frac{\kappa \tilde{q}}{2} s^\rho \epsilon _{\sigma \rho \pi }-\frac{1}{2}(u^\delta \epsilon _{\sigma \rho \pi }\nabla ^\rho \nabla _\delta f_R), \end{aligned}$$
(35)

with \(\psi _{\sigma \pi }\) and \(\varphi _{\sigma \pi }\) being equal to

$$\begin{aligned} \psi _{\sigma \pi }&=\frac{h_{\sigma \theta }\nabla ^\theta \nabla _\pi f_R }{4}-\frac{(f-Rf_R)h_{\sigma \pi }}{6}\\&\quad +\frac{u_\pi u^\epsilon \nabla _\sigma \nabla _\epsilon f_R}{4}-\frac{u^\rho u_\pi \nabla _\sigma \nabla _\rho f_R}{4}\\&-\frac{u_\sigma u_\theta \nabla ^\theta \nabla _\pi f_R}{4}-\frac{(\Box f_R)g_{\sigma \pi }}{4}\\&\quad +\frac{(\Box f_R)u_\sigma u_\pi }{2}+\frac{\nabla _\sigma \nabla _\pi f_R}{4}\\ {}&+\frac{\delta _{\sigma \pi } u_\theta u^\epsilon \nabla ^\theta \nabla _\epsilon f_R}{8}+\frac{\delta _{\sigma \pi } u_\theta u^\rho \nabla ^\theta \nabla _\rho f_R}{8}\\&\quad -\frac{\delta _{\sigma \pi } h^\rho _\theta \nabla ^\theta \nabla _\rho f_R}{8}-\frac{\delta _{\sigma \pi } h^\epsilon _\theta \nabla ^\theta \nabla _\epsilon f_R}{8},\\ \varphi _{\sigma \pi }&=\frac{1}{2}\left[ \nabla _\sigma \nabla _\pi f_R-u_\pi u^\delta \nabla _\sigma \nabla _\delta f_R\right. \\&\quad \left. -u_\sigma u_\theta \nabla ^\theta \nabla _\pi f_R +g_{\sigma \pi }u_\theta u^\delta \nabla ^\theta \nabla _\delta f_R\right] \\&\quad +\left( \frac{f-Rf_R}{6}\right) h_{\sigma \pi }. \end{aligned}$$

Utilizing the splitting of the Bel tensor [39], we acquire Bel superenergy \(\bar{W}\) and super-Poynting vector \(\bar{P}_\sigma \). \(\bar{W}\), unlike its name, does not express energy rather it is sometimes used to explain the existence of vorticity in radiative and stationary metrics. Super-Poynting vector \(\bar{P}_\sigma \) follows the same concept as that of the Poynting vector in electromagnetism. It indicates the flux of the superenergy and assists in comprehending gravitational radiation. The clear physical meaning of both of these quantities is still a mystery. However, the expressions for both are given as follows

$$\begin{aligned} \bar{W}&=\frac{1}{2}\left( X_{\sigma \pi }X^{\sigma \pi }+Y_{\sigma \pi }Y^{\sigma \pi }\right) +Z_{\sigma \pi }Z^{\sigma \pi },\\ \bar{P}_\sigma&=\epsilon _{\sigma \pi \gamma }\left( Y^\gamma _\delta Z^{\pi \delta }-X^\gamma _\delta Z^{\delta \pi }\right) . \end{aligned}$$

Using the values of \(X_{\sigma \pi },~Y_{\sigma \pi }\) and \(Z_{\sigma \pi }\), we obtain the following outcome

$$\begin{aligned} \bar{W}&=\frac{5\kappa ^2 \tilde{\rho }^2}{24}+\frac{\kappa ^2 \Pi ^2}{6}+\frac{2E^2}{3}+\frac{\kappa ^2 \tilde{\rho }\check{P}}{4}\nonumber \\ {}&\quad +\frac{3\kappa ^2 \check{P}^2}{8}+\frac{H_1+H_2}{2}+\frac{\kappa ^2 \tilde{q}^2}{2}+\psi _1, \end{aligned}$$
(36)
$$\begin{aligned} \bar{P}_\sigma&=\frac{\kappa ^2 \tilde{q}}{2}\left( \tilde{\rho }+\tilde{P}_r\right) s_\sigma +Y^\gamma _\delta \vartheta ^\delta _{\sigma \gamma }-X^\gamma _\delta \vartheta ^\delta _{\sigma \gamma } \end{aligned}$$
(37)

with \(H_1,H_2\) and \(\vartheta ^\delta _{\sigma \gamma }\) being equal to

$$\begin{aligned} H_1&=\psi _{\sigma \pi }\times \left( \frac{\kappa \tilde{\rho }}{3}h^{\sigma \pi }+\frac{\kappa \Pi ^{\sigma \pi }}{2}-E^{\sigma \pi } +\psi ^{\sigma \pi }\right) ,\\ H_2&=\varphi _{\sigma \pi } \times \left( \frac{\kappa \tilde{\rho }}{3}h^{\sigma \pi }+\frac{\kappa \Pi ^{\sigma \pi }}{2}-E^{\sigma \pi }+\varphi ^{\sigma \pi }\right) ,\\ \vartheta ^\delta _{\sigma \gamma }&=-\frac{1}{2}\epsilon _{\sigma \pi \gamma }\epsilon ^{\pi \delta }_\rho u^\theta \nabla ^\rho \nabla _\theta f_R. \end{aligned}$$

An expression for superenergy found using the Bel–Robinson tensor given as \(W=E^{\sigma \pi }E_{\sigma \pi }\) becomes \(W=\frac{2E^2}{3}\). It can be seen that W signifies the effects of tidal forces in that specific stellar equation. We evaluated an expression for \(\bar{W}-W\) to link the fluid parameters with the superenergy excluding the effects of tidal forces. Combining this result along with Eq. (36), we attain

$$\begin{aligned} \bar{W}-W&=\frac{5\kappa ^2 \tilde{\rho }^2}{24}+\frac{\kappa ^2 \Pi ^2}{6}+\frac{\kappa ^2 \tilde{\rho }\check{P}}{4}\\&\quad +\frac{3\kappa ^2 \check{P}^2}{8}+\frac{H_1+H_2}{2}+\frac{\kappa ^2 \tilde{q}^2}{2}+\psi _1,\\ \vartheta ^\delta _{\sigma \gamma }&=-\frac{1}{2}\epsilon _{\sigma \pi \gamma }\epsilon ^{\pi \delta }_\rho u^\theta \nabla ^\rho \nabla _\theta f_R. \end{aligned}$$

3.1 Five scalar functions

This section focusses on the derivation of five scalar entities termed as structure scalars which are further used in the delineation of dynamics of the astrophysical objects. First, we split \(X_{\sigma \pi }\) and \(Y_{\sigma \pi }\) into their trace and trace-free parts as follows

$$\begin{aligned} X_{\sigma \pi }=\frac{1}{3}Tr X h_{\sigma \pi }+X_{<\sigma \pi >}, \end{aligned}$$

with \(Tr X=X^\alpha _\alpha \) so that

$$\begin{aligned} X_{\langle \sigma \pi \rangle }=h^\mu _\sigma h^\nu _\pi \left( X_{\mu \nu }-\frac{1}{3} Tr X h_{\mu \nu }\right) , \end{aligned}$$

with the trace part TrX of \(X_{\sigma \pi }\) being equal to

$$\begin{aligned} X_T&=\left[ \kappa \tilde{\rho }-\left( \frac{f-Rf_R}{2}\right) \right. \nonumber \\&\quad \left. -\frac{h_\pi ^\epsilon \nabla ^\pi \nabla _\epsilon f_R}{4}-\frac{h^\rho _\pi \nabla ^\pi \nabla _\rho f_R}{4}\right. \nonumber \\&\quad \left. -\frac{h^\epsilon _\theta \nabla ^\theta \nabla _\epsilon f_R}{4}-\frac{h^\rho _\theta \nabla ^\theta \nabla _\rho f_R}{4}\right] . \end{aligned}$$
(38)

It is worth observing that

$$\begin{aligned} X_{\langle \sigma \pi \rangle }=X_{TF}\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) , \end{aligned}$$

where \(X_{TF}\) comes out to be as follows

$$\begin{aligned} X_{TF}\equiv \left( \frac{\kappa \Pi }{2}-E+\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) } \right) , \end{aligned}$$
(39)

Following the same pattern, we acquire the trace-free part \(Y_T\) of \(Y_{\sigma \pi }\)

$$\begin{aligned} Y_T&=\left[ \frac{\kappa }{2} (\tilde{\rho }+3\tilde{P}_r-2\Pi )\right. \nonumber \\&\quad \left. +\frac{f-Rf_R}{2}+\frac{\Box f_R}{2}-\frac{u^\sigma u^\delta \nabla _\sigma \nabla _\delta f_R}{2}\right. \nonumber \\&\quad \left. -\frac{u^\pi u_\theta \nabla ^\theta \nabla _\pi f_R}{2}+2u_\theta u^\delta \nabla ^\theta \nabla _\delta f_R\right] . \end{aligned}$$
(40)

We can express

$$\begin{aligned} Y_{\langle \sigma \pi \rangle }=Y_{TF}\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) , \end{aligned}$$

with the value of \(Y_{TF}\) as follows

$$\begin{aligned} Y_{TF}\equiv \left( \frac{\kappa \Pi }{2}+E+\frac{\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) } \right) . \end{aligned}$$
(41)

The last scalar is found by using the tensor \(Z_{\sigma \pi }\) as

$$\begin{aligned} Z=\sqrt{Z_{\sigma \pi }Z^{\sigma \pi }}=\sqrt{\left( \frac{\kappa ^2 \tilde{q}^2}{2}+\psi _1\right) }. \end{aligned}$$
(42)

The scalar functions \(X_{TF}\) and \(Y_{TF}\) are found to correspond with the pressure anisotropy along with extra curvature terms of f(R) gravity theory as follows

$$\begin{aligned} \kappa \Pi =X_{TF}+Y_{TF}-\frac{\xi _{\sigma \pi }+\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }. \end{aligned}$$
(43)

Using these five scalar functions, we rewrite structure and evolution equations so that one can easily observe which physical characteristic of the fluid as described by each equation. Thus, Eqs. (24)–(31) can be restated as

$$\begin{aligned}&\frac{\kappa \tilde{\rho }^*}{2}+\frac{1}{3}\left[ X_{T}+X_{TF}+Y_{T}+Y_{TF}-\chi _4\right. \nonumber \\&\quad \left. -\chi _3-\frac{\xi _{\sigma \pi }+\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right] \theta = \frac{1}{3}\left( \theta +\frac{\sigma }{2}\right) \nonumber \\&\quad \times \left( X_{TF}+Y_{TF}+\frac{\xi _{\sigma \pi }+\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right) \nonumber \\&\quad -\frac{\sqrt{2}}{2}\left( \sqrt{Z^2-\psi _1}\right) ^\dagger -\sqrt{2}\left( \sqrt{Z^2-\psi _1}\right) \nonumber \\&\quad \times \left[ a+\frac{s^1}{r}\right] -\frac{\kappa \chi _1}{2}, \end{aligned}$$
(44)
$$\begin{aligned}&\frac{\kappa \tilde{P}_r^\dag }{2}+\frac{a}{3}\left[ X_{T}+X_{TF}+Y_{T}+Y_{TF}-\chi _4-\chi _3\right. \nonumber \\&\qquad \left. -\frac{\xi _{\sigma \pi }+\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right] + \frac{s^1}{r}\nonumber \\&\qquad \times \left( X_{TF}+Y_{TF}+\frac{\xi _{\sigma \pi }+\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right) \nonumber \\&\quad =\frac{\sqrt{2}}{3}\left( \sqrt{Z^2-\psi _1}\right) \left( \frac{\sigma }{2}-2\theta \right) -\frac{\sqrt{2}}{2} \nonumber \\&\qquad \times \left( \sqrt{Z^2-\psi _1}\right) ^*-\frac{\kappa \chi _2}{2}, \end{aligned}$$
(45)
$$\begin{aligned}&\theta ^*+\frac{\theta ^2}{3}+\frac{\sigma ^2}{6}-a^\dagger -a^2-\frac{2as^1}{r}=-Y_T, \end{aligned}$$
(46)
$$\begin{aligned}&\left( \frac{\sigma }{2}+\theta \right) ^\dagger =-\frac{3\sigma s^1}{2r}+\frac{3\sqrt{2}Z}{2}, \end{aligned}$$
(47)
$$\begin{aligned}&a^\dagger +a^2+\frac{\sigma ^*}{2}+\frac{\theta \sigma }{3}-\frac{as^1}{r}-\frac{\sigma ^2}{12}=-Y_{TF}, \end{aligned}$$
(48)
$$\begin{aligned}&\frac{1}{3f_R}\left[ Y_{T}-Y_{TF}-2X_{TF}+X_{T}-\chi _4-\chi _3\right. \nonumber \\&\qquad \left. -\frac{2\xi _{\sigma \pi }+3\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right] \left( \frac{\sigma }{2}+\theta \right) \nonumber \\&\qquad +\left[ \frac{X_{T}}{2}-X_{TF}-\frac{\chi _4}{2}+\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right] ^*\nonumber \\&\quad =-\frac{3s^1}{2r}(\sqrt{2Z^2-\psi _1}), \end{aligned}$$
(49)
$$\begin{aligned}&\left( \frac{\kappa \tilde{\rho }}{2}+\frac{\chi _4}{2}-X_{TF}\right) ^\dagger =\frac{3s^1 }{r}X_{TF}+\frac{\sqrt{2}Z}{2}\left( \frac{\sigma }{2}+\theta \right) , \end{aligned}$$
(50)
$$\begin{aligned}&\frac{3m}{r^3}=\frac{1}{f_R}\left[ \frac{X_{T}}{2}-X_{TF}-\frac{\chi _4}{2} -\frac{\tau ^0_0+\tau ^1_1}{2}+\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right] . \end{aligned}$$
(51)

We can also express Bel superenergy and Super-Poynting vector in terms of these five scalar functions. Thus, Eq. (36) and (37) can be rewritten as

$$\begin{aligned} \bar{W}&=\frac{1}{6}\left[ X_{T}^2+Y_{T}^2-D_3\right] \\&\quad +\frac{1}{3}\left[ X_{TF}^2+Y_{TF}^2+D_4\right] +Z^2-\psi _1,\\ \bar{P}_\sigma&=\frac{\sqrt{2Z^2-\psi _1}}{3}\left[ X_{T}+Y_{T}+D_1+D_2+X_{TF}\right. \\&\quad \left. +Y_{TF}+\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right] s_\sigma , \end{aligned}$$

where we have

$$\begin{aligned} D_1&=-\left( \frac{f-Rf_R}{2}\right) -\frac{h_\pi ^\epsilon \nabla ^\pi \nabla _\epsilon f_R}{4}-\frac{h^\rho _\pi \nabla ^\pi \nabla _\rho f_R}{4}\\&\quad -\frac{h^\epsilon _\theta \nabla ^\theta \nabla _\epsilon f_R}{4}-\frac{h^\rho _\theta \nabla ^\theta \nabla _\rho f_R}{4},\\ D_2&=\frac{f-Rf_R}{2}+\frac{\Box f_R}{2}-\frac{u^\sigma u^\delta \nabla _\sigma \nabla _\delta f_R}{2}\\&\quad -\frac{u^\pi u_\theta \nabla ^\theta \nabla _\pi f_R}{2}+2u_\theta u^\delta \nabla ^\theta \nabla _\delta f_R,\\ D_3&=D_1^2+\kappa \tilde{\rho }D_1+D_2^2+\kappa \rho D_2+3\kappa \check{P}D_2,\\ D_4&=2E\left( \frac{\xi _{\sigma \pi }+\zeta _{\sigma \pi }}{s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}}\right) \\&\quad +\kappa \Pi \left( \frac{\zeta _{\sigma \pi }-\xi _{\sigma \pi }}{s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}}\right) +\left( \frac{\xi _{\sigma \pi }}{s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}}\right) ^2\\&\quad +\left( \frac{\zeta _{\sigma \pi }}{s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}}\right) ^2 \end{aligned}$$

3.2 Physical interpretation of structure scalars

In this section, we interpret five scalar functions by comprehending the physical processes connected with these scalars. It is evident that \(X_T\) deals with the energy density of the fluid under consideration while Z provides information about the dissipative flux of the fluid. Both these scalars incorporate extra curvature terms delineating the fact that we are working in f(R) theory of gravity. To learn about the physical interpretation of \(Y_T\) and \(Y_{TF}\), we get back to the expression for Tolman mass to get the following result

$$\begin{aligned} m_T&=(m_T)_\Sigma \left( \frac{r}{r_\Sigma }\right) ^3\\&\quad -r^3\int _r^{r_\Sigma }\frac{r^3e^{(\nu -\lambda )/2}}{2r}\left( \nu 'f_R'\right) dr\\&\quad -r^3\int _r^{r_\Sigma }\frac{f_Re^{(\lambda -\nu )/2}}{2r}\ddot{\lambda }dr \\&\quad + r^3\int _r^{r_\Sigma }\frac{e^{(\nu +\lambda )/2}}{r}\left[ \left( Y_{TF}-\frac{\chi _2}{\left( s_\alpha s_\beta +\frac{h_{\alpha \beta }}{3}\right) }\right. \right. \\&\quad \left. \left. -\frac{\kappa }{2}\left( \left. T^1_1\right. ^{(D)}-2\left. T^2_2\right. ^{(D)} \right) \right) +Ef_R\right] dr. \end{aligned}$$

We observe that the Tolman mass expression for a system in equilibrium or quasi-equilibrium becomes

$$\begin{aligned} m_T=\int _0^rr^2 e^{(\lambda +\nu )/2}(Y_T-\chi _3)dr. \end{aligned}$$

4 Static spheres anisotropic in pressure

The center of attention of this section is to express few line elements corresponding to static spheres as an application of our work. Three different formalisms are used to arrive at such spherically symmetric line elements bearing in mind the f(R) theory of gravity. Such spheres are useful in comprehending the evolution of static astrophysical objects. Moreover, we utilize structure scalars in the stellar equations of such static objects to comprehend the physical processes explained by those equations.

4.1 First alternative

Utilizing Eqs. (13) and (51), we acquire the following outcome

$$\begin{aligned} e^{-\lambda }= & {} 1-\frac{2r^2}{3f_R}\left[ \frac{X_{T}}{2}-X_{TF}-\frac{\chi _4}{2}\right. \nonumber \\&\left. -\frac{\tau ^0_0+\tau ^1_1}{2}+\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right] . \end{aligned}$$
(52)

Equations (46) and (48) can be restated for static case as

$$\begin{aligned} a=\frac{r}{3s^1}(Y_{TF}+Y_T). \end{aligned}$$
(53)

We also have another relation as follows

$$\begin{aligned} a=\frac{e^{-\lambda /2}\nu '}{2};\quad \quad \quad s^1=e^{-\lambda /2}. \end{aligned}$$
(54)

Putting Eq. (54) into Eq. (53) and performing integration w.r.t the radial coordinate, we attain

$$\begin{aligned} e^\nu =Ce^{\int \frac{2r}{3}(Y_{TF}+Y_T)\left[ 1-\frac{2r^2}{3f_R}\left( \frac{X_{T}}{2}-X_{TF}-\frac{\chi _4}{2} -\frac{\tau ^0_0+\tau ^1_1}{2}+\frac{\xi _{\sigma \pi }}{s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}}\right) \right] dr}, \end{aligned}$$
(55)

with C indicating the integration constant. Making use of Eq. (52) and (55), the static spherically symmetric line element in this case becomes

$$\begin{aligned} ds^2&=Ce^{\int \frac{2r}{3}(Y_{TF}+Y_T)\left[ 1-\frac{2r^2}{3f_R}\left( \frac{X_{T}}{2}-X_{TF}-\frac{\chi _4}{2} -\frac{\tau ^0_0+\tau ^1_1}{2}+\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right) \right] dr} dt^2\nonumber \\&\quad -\left[ 1-\frac{2r^2}{3f_R}\left( \frac{X_{T}}{2}-X_{TF}-\frac{\chi _4}{2} -\frac{\tau ^0_0+\tau ^1_1}{2}\right. \right. \nonumber \\&\quad \left. \left. +\frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right) \right] ^{-1}dr^2-r^2 sin^2 \theta d\phi ^2. \end{aligned}$$
(56)

4.2 Second alternative

Another approach to arrive at the line element is by making use of Eqs. (15), (38) and (51) as follows

$$\begin{aligned} m(r)=\frac{r^3}{3}\left[ \frac{m'}{r^2}-\frac{X_{TF}}{f_R}+\frac{\chi _4+\tau ^0_0}{2f_R}+M_1\right] , \end{aligned}$$

with \(M_1\) being equal to

$$\begin{aligned} M_1=\frac{1}{f_R}\left[ \frac{\xi _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }-\frac{\chi _4}{2} -\frac{\tau ^0_0+\tau ^1_1}{2}\right] . \end{aligned}$$

Upon integration, the following outcome is attained

$$\begin{aligned} m(r)=r^3\int \left( \frac{X_{TF}}{rf_R}-\frac{\chi _4+\tau ^0_0}{2rf_R}-\frac{M_1}{r}\right) dr+c_1. \end{aligned}$$

Utilizing Eq. (13) here, we deduce

$$\begin{aligned} e^{-\lambda }=1-2r^2\left[ \int \left( \frac{X_{TF}}{rf_R}-\frac{\chi _4+\tau ^0_0}{2rf_R}-\frac{M_1}{r}\right) dr+c_1\right] . \end{aligned}$$

Keeping in view the static case, the modified field equation (3) can be restated as

$$\begin{aligned} \frac{\kappa P_r}{2}=\frac{f_R}{2}\left( \frac{e^{-\lambda }-1}{r^2}\right) +\frac{\nu 'e^{-\lambda }f_R}{2r}+\frac{\tau ^1_1}{2}. \end{aligned}$$

Next, we utilize Eq. (13), (40), (41), (53) and (54)

$$\begin{aligned} \frac{\kappa P_r}{2}+\frac{mf_R}{r^3}-\frac{\tau ^1_1}{2}=\frac{\nu 'e^{-\lambda }f_R}{2r}=Y_h, \end{aligned}$$
(57)

with \(Y_h\) being

$$\begin{aligned} Y_h=\frac{f_R}{3}\left( Y_T+Y_{TF}\right) +\frac{\tau ^1_1}{2}. \end{aligned}$$

Performing integration of Eq. (57), we obtain

$$\begin{aligned} e^\nu =C_2e^{\int \frac{2r}{f_R}Y_h\left[ 1-2r^2\left( \int \left( \frac{X_{TF}}{rf_R}-\frac{\chi _4+ \tau ^0_0}{2rf_R}-\frac{M_1}{r}\right) dr+c_1\right) \right] ^{-1}dr}, \end{aligned}$$
(58)

with \(C_2\) indicating the integration constant. The line element takes the form

$$\begin{aligned} ds^2&=C_2e^{\int \frac{2r}{f_R}Y_h\left[ 1-2r^2\left( \frac{X_{TF}}{rf_R}-\frac{\chi _4+\tau ^0_0}{2rf_R} -\frac{M_1}{r}\right) dr+c_1\right] ^{-1}dr}dt^2\nonumber \\&-\left[ 1-2r^2 \int \left( \frac{X_{TF}}{rf_R} -\frac{\chi _4+\tau ^0_0}{2rf_R}-\frac{M_1}{r}\right) dr+c_1\right] ^{-1}dr^2\nonumber \\&-r^2d\theta ^2-r^2 sin^2 \theta d\phi ^2. \end{aligned}$$
(59)

Both of these alternative approaches produced line elements by utilizing only a combination of four structure scalars. In the next subsection, we use another formalism developed in [40] to deduce third line element that utilizes only two scalars to produce the required line element.

4.3 Third alternative

This formalism utilizes the scalars \(X_{TF}\) and \(Y_{TF}\) to obtain the spherically symmetric line element. From the electric part of Weyl tensor, we may write

$$\begin{aligned} E=-\frac{e^{-\lambda }}{2}\left[ \frac{\nu ''}{2}+\left( \frac{\nu '}{2}\right) ^2+ \frac{\nu '}{2}\left( -\frac{\lambda '}{2}-\frac{1}{r}\right) +\frac{\lambda '}{2r}+ \frac{1-e^\lambda }{r^2}\right] . \end{aligned}$$
(60)

Now, writing new variables y and u as given in [40]

$$\begin{aligned} y=e^{-\lambda }; \quad \quad \quad \frac{\nu '}{2}=\frac{u'}{u}. \end{aligned}$$

Equation (60) can then be restated as

$$\begin{aligned} y'+2y\left( \frac{u''-\frac{u'}{r}+\frac{u}{r^2}}{u'-\frac{u}{r}}\right) =\frac{2u(1-2r^2E)}{r^2\left( u'-\frac{u}{r}\right) }. \end{aligned}$$

Upon integrating w.r.t. r, we get

$$\begin{aligned} y=e^{-\int k(r)dr}\left( \int e^{-\int k(r)dr}f(r)dr+c_1\right) , \end{aligned}$$

with k(r) and f(r) as follows

$$\begin{aligned} k(r)=2 \frac{d}{dr}\left[ ln\left( u'-\frac{u}{r}\right) \right] , \quad \quad f(r)=\frac{2u(1-2r^2E)}{r^2\left( u'-\frac{u}{r}\right) }. \end{aligned}$$

The alphabet \(c_1\) indicates the constant of integration procured using the matching conditions. Utilizing original variables, we attain

$$\begin{aligned}&\frac{\nu '}{2}-\frac{1}{r}=\frac{e^{\lambda /2}}{r}\nonumber \\&\quad \times \sqrt{(1-2Er^2)+c_1r^2e^{-\nu }+ r^2e^{-\nu }\int \frac{e^\nu }{r^2}(2r^2E)'dr}.\nonumber \\ \end{aligned}$$
(61)

Making use of a new variable z introduced in [40], we write

$$\begin{aligned} e^\nu =\frac{e^{2\int zdr}}{r^2}, \end{aligned}$$
(62)

with the help of which following outcome is established

$$\begin{aligned} z(r)=\frac{\nu '}{2}+\frac{1}{r}. \end{aligned}$$
(63)

Utilizing Eqs. (62) and (63) to create a link between z and E as follows

$$\begin{aligned}&z(r)=\frac{2}{r}+\frac{e^{\lambda /2}}{r}\sqrt{(1-2Er^2)+c_1r^4e^{-\int 2z(r)dr}+ r^4e^{-\int 2z(r)dr}\int \frac{e^{-\int 2z(r)dr}}{r^4}(2r^2E)'dr}. \end{aligned}$$
(64)

Modified field equations provide us with the following outcome

$$\begin{aligned} \kappa \Pi&=f_R\left[ e^{-\lambda }\left( -\frac{\nu ''}{2}-\left( \frac{\nu '}{2}\right) ^2+ \frac{\nu '}{2r}+\frac{1}{r^2}\right) \right. \\&\quad \left. +\frac{\lambda 'e^{-\lambda }}{2}\left( \frac{\nu '}{2}+ \frac{1}{r}\right) -\frac{1}{r^2}\right] +(\tau ^2_2-\tau ^1_1). \end{aligned}$$

where

$$\begin{aligned} \tau ^2_2= & {} \ddot{f_R}e^{-\nu }-\frac{f_R'e^{-\lambda }}{r}+\frac{\dot{\nu }\dot{f_R} e^{-\nu }}{2}-\frac{\nu ' f_R'e^{-\lambda }}{2}-f_R'' e^{-\lambda }\\&+\frac{\dot{\lambda } \dot{f_R} e^{-\nu }}{2}+\frac{\lambda 'f_R'e^{-\lambda }}{2}-\frac{(f-Rf_R)}{2}. \end{aligned}$$

Using z and y, we can write

$$\begin{aligned}&y'+y\left[ \frac{2z'}{z}+2z-\frac{6}{r}+\frac{4}{r^2z}\right] \\&\quad =-\frac{2}{z}\left[ \frac{\kappa \Pi }{f_R} +\frac{1}{r^2}-\frac{(\tau ^2_2-\tau ^1_1)}{f_R}\right] . \end{aligned}$$

Integrating using the radial coordinate r lead us to the following outcome

$$\begin{aligned} e^{\lambda (r)}=\frac{z^2 e^{\int \left( 2z+\frac{4}{zr^2}\right) dr}}{r^6\left[ -2\int \frac{z}{r^8}\left( \frac{\kappa \Pi r^2}{f_R}+1-\frac{(\tau ^2_2-\tau ^1_1)r^2}{f_R}\right) e^{\int \left( 2z+\frac{4}{zr^2}\right) dr}dr+C\right] }, \end{aligned}$$
(65)

with C being the integration constant. Utilizing the scalar functions \(X_{TF}\) and \(Y_{TF}\), Eqs. (64) and (65) can be used to establish the following result

$$\begin{aligned} z(r)&=\frac{2}{r}+\frac{e^{\lambda /2}}{r}\left[ 1-r^2\left( Y_{TF}-X_{TF}+\frac{\xi _{\sigma \pi }-\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) }\right) \right. \nonumber \\&\quad + r^4e^{-\int 2zdr}\left( c_1+\int \frac{e^{\int 2zdr}}{r^4}\right. \nonumber \\&\left. \left. \times \left( r^2(Y_{TF}-X_{TF}+\frac{\xi _{\sigma \pi }-\zeta _{\sigma \pi }}{\left( s_\sigma s_\pi +\frac{h_{\sigma \pi }}{3}\right) })\right) 'dr\right) \right] ^{1/2}, \end{aligned}$$
(66)

and

$$\begin{aligned} e^{\lambda (r)}=\frac{z^2 e^{\int \left( 2z+\frac{4}{zr^2}\right) dr}}{r^6\left[ -2\int \frac{z}{r^8} \left( \frac{r^2}{f_R}+r^2(Y_{TF}+X_{TF})\right) e^{\int \left( 2z+\frac{4}{zr^2}\right) dr}dr+C\right] }. \end{aligned}$$

If we take into account the conformally flat fluids anisotropic in pressure, the scalar \(Y_{TF}\) becomes equal to \(X_{TF}\) plus extra curvature terms.

5 Conclusion

In this manuscript, we have looked into the f(R) corrections on the dynamical attributes of stellar configurations. To carry out this study, we have taken into account the anisotropic spherically symmetric dissipative geometry coupled with radiating matter constituents. We evaluated certain kinematical quantities including shear tensor and expansion scalar to comprehend the basic fluid properties. We then manipulated the modified field equations keeping in view the f(R) theory of gravitation. Matching conditions are attained in order to smoothly match the interior and exterior geometry. The tensor quantity responsible for the production of tidal forces, i.e., the Weyl tensor along with the Riemann tensor is also evaluated which are later on linked with matter distribution and metric variables with the addition of dark source terms arising due to f(R) theory of gravitation. Moreover, the expression for Misner–Sharp mass function is calculated to learn about the quantity of matter lying inside spherically symmetric geometrical object. Interesting relationships between Misner–Sharp mass, metric variables and the Weyl scalar equation are also established because of their extraordinary significance in the modeling of stellar configurations. Few structure and evolution stellar equations including Bianchi identities, Ricci identities and Raychaudhuri equation for the study of dynamics of stellar objects are also computed.

Further, we split the expression for the Riemann curvature tensor by utilizing the orthogonal decomposition formalism introduced by Herrera et al. [40]. We utilized this formalism on our stellar model in the background of f(R) theory of gravitation. We found expressions for three tensors namely, \(X_{\sigma \pi },~Y_{\sigma \pi }\) and \(Z_{\sigma \pi }\). We then manipulated the trace and trace-free parts of these three tensors due to their peculiar significance in the analysis of gravitational collapse and stellar evolution. Such trace and trace-free parts are scalar functions termed as f(R) structure scalars. We express the previously computed structure and evolution equations in terms of these f(R) structure scalars and the resulting equations are found to be extremely useful in the study of Penrose–Hawking singularity theorems. The part and parcel of such theorems is the Raychaudhuri equation which signifies the dynamics of flow in spacetime. Utilizing the f(R) structure scalars in Raychaudhuri equation and then using it in those theorems allows us to comprehend the occurrence of singularities in spacetime. That’s how the f(R) structure scalars are beneficial in the interpretation of Penrose–Hawking singularity theorems.

It is worth noticing that the scalar function \(X_T\) deals with the energy density of the fluid under consideration with some extra curvature terms added in its expression due to the effect of f(R) theory of gravitation. Theses extra curvature ingredients symbolize the presence of unknown form of matter or energy, i.e., DM or DE. The expression for \(X_{TF}\) and \(Y_{TF}\) comprises pressure anisotropy in addition to electric part of Weyl tensor and extra curvature terms. \(Y_T\) provides information about the pressure anisotropy, energy density and few terms expressing the effects of DM or DE.

As a practical example of our work, we may contemplate a toy model \(f(R)= R +\varepsilon R^2\) known as Starobinsky model. Since different models characterize different evolutionary phases of our Universe, Starobinsky model signifies the inflationary era. The extra curvature quantities of this model may unravel the mystery of influence of DM or DE on the kinematics of our Universe. Utilization of structure scalars and this value of f(R) in the stellar equations obtained in this manuscript would describe the effects of modified gravity on the matter distribution under consideration.