Raychaudhuri-based reconstruction of anisotropic Einstein-Maxwell equation in 1+3 covariant formalism of $f(R)$-gravity

Recently, a new strategy to the reconstruction of $f(R)$-gravity models based on the Raychaudhuri equation has been suggested by Choudhury et al. In this paper, utilizing this method, the reconstruction of anisotropic Einstein-Maxwell equation in the $1+3$ covariant formalism of $f(R)$-gravity is investigated in four modes: $i.$ Reconstruction from a negative constant deceleration parameter refereeing to an ever-accelerating universe; $ii.$ Reconstruction from a constant jerk parameter $j=1$ which recovers celebrated $\Lambda \text{CDM}$ mode of evolution; $iii.$ Reconstruction from a variable jerk parameter $j=Q(t)$; and $iv.$ Reconstruction from a slowly varying jerk parameter. Furthermore, two suggestions for enhancing the method are proposed.

It is argued that an unknown hiddenly characteristics sort of energy with large negative pressure is responsible for this accelerating cosmic expansion [8]. This mysterious candidate being incompatible with strong energy condition is dubbed as dark energy. This cosmic behavior can also be explained either by modifying matter part ( f (R) gravity, f (T ) gravity, f (T ) gravity with unusual term [9], scalar-tensor theories [10], etc.) or by modifying geometric part (chaplygin gas [11], quintessence [12], phantom [13], quintom [14], etc.) of the Einstein-Hilbert action. This new set of gravity theories passes several solar system and astrophysical tests successfully [15,16].
The simple modification of Einstein's theory of gravity namely f (R) -gravity, as a source of acceleration, was proposed by Capozziello et al. [17] and Carroll et al. [18]. The f (R) model gives sufficient generality to encapsulate some of the basic characteristics of higher-order gravity and yet are rather simple to handle. The modified f (R) -gravity can elucidate the cosmic acceleration without introducing the dark energy component [19][20][21]. In addition, it has been demonstrated that the modified f (R) -gravity can be derived from string/M-theory [22].
One of the unsolved problems in cosmology is the cosmological magnetic field emerging at large scale of the universe observationally [23]. In order to disentangle the origin of the primordial cosmological magnetic field, there are many theoretical explanations, for example, it has been created from the Big Bang like all matters populating the universe [24]. Besides this, such fields might play some roles on the cosmic microwave background radiation. For these reasons, in this paper, the primordial magnetic fields are included in the energy-momentum tensor of the Einstein field equation directly. It is worthwhile mentioning that the cosmological magnetic fields will naturally appear in the universe when the anisotropic cosmological models are taken into account. Hence we would like to consider the problem in an anisotropic background.
In this paper, from the perspective of this new method, we investigate the reconstruction of anisotropic Einstein-Maxwell equation in 1+3 covariant formalism of f (R) -gravity in four modes of evolution: i. Reconstruction from a negative constant deceleration parameter refereeing to an ever-accelerating universe; ii. Reconstruction from a constant jerk parameter j = 1 which recovers celebrated ΛCDM mode of evolution; iii. Reconstruction from a time-variable jerk parameter j = Q(t) ; and iv. Reconstruction from a slowly varying jerk parameter. The two last jerk types for f (T ) -gravity has recently been studied by Chakrabarti et al. [39]. Finally, some suggestions for enhancing the method are proposed.

The model and basic equations
In this section, the evolution equations of the f (R) -gravity in orthogonally spatially homogeneous 1+3 covariant approach have been set up.
For a given fluid four-velocity vector field u µ , the projection tensor h µν = g µν + u µ u ν projects into the instantaneous rest-space of a comoving observer who, in this paper, is characterized by u µ = (1, 0, 0, 0) . Indeed, the four-velocity u µ is orthogonal to the induce metric h µν (i.e. h µν u µ = 0 ).
Introducing the vorticity tensor ω µν ( ω µν = ω [µν] , ω µν u ν = 0 ), the symmetric shear tensor σ µν ( σ µν = σ (µν) , σ µν u ν = 0 , σ α α = 0 ), and the volume expansion Θ = ∇ α u α the first covariant derivative of the four-velocity can therefore be decomposed as whereu µ , which is defined asu µ = u ν ∇ ν u µ , is the acceleration vector. The last term in (1) is indeed the following difference where Θ µν are the components of the volume expansion tensor of the fluid (or the extrinsic curvature) whose its trace (i.e. Θ ≡ Θ µν h µν ) is the rate of the volume expansion parameter namely Hubble parameter. Relatively to u µ , the energy-momentum tensor can be decomposed in the form: where ρ is the energy density, p is the isotopic pressure, q µ is the energy flux ( q α u α = 0 ), and π µν is the symmetric trace-free anisotropic stress pressure ( π µν = π νµ , π α α = 0 , π µν u µ = 0 ). We start with the gravitational action of f (R) -gravity of the form: where g is the determinant of metric, f (R) is a function of the Ricci scalar R , and L M stands for the matter fields Lagrangian density. Varying this action with respect to metric tensor g µν , the Einstein field equations are obtained as follows 1 : where f (R) = df /dR , and all the subscripts and superscripts run from zero to three (i.e. 0 , 1 , 2 , and 3 ). Assuming the total energy-momentum tensor, T M µν , consists of an electromagnetic filed, T em µν , and a perfect fluid, T pf µν , as two main non-interacting parts, it can therefore be written as and in which F µν is the field strength. For a given electric field, E µ , and magnetic field, B µ , the field strength F µν is defined as in which η µναβ is an antisymmetric permutation tensor of space-time with η 0123 = √ −g . The energy-momentum tensor of Maxwell field can be recast in the form where ρ em and p em are the energy density and the isotropic pressure of the electromagnetic field, respectively, and they are given by and the anisotropic stress is In the present paper, we prefer to work with a pure magnetic case (i.e. E = 0 and B = 0 ). Let us consider the problem in the anisotropic background of the form This line element is known as Locally Rotationally Symmetric Bianchi type-I (LRS B-I). Pursuant to this background geometry, the magnetic fields may have component as B µ = (0; B(t), 0, 0) . Defining S µν = ∇ µ ∇ ν f and using equation (5), the Ricci tensor takes the following form Now, utilizing this equation, the Ricci tensor R µν can be split into the following forms: where p = p tot. = p m + p em . In analogous with the Ricci tensor, for the S µν one has the following relations: Therefore, it can be demonstrated that the Raychaudhuri equation in the 1+3 covariant formalism of f (R) -gravity of Bianchi type-I is obtained asΘ where ρ = ρ tot. = ρ m + ρ em . Restricting the magnetic field to be aligned along the shear eigenvector, the shear tensor would also diagonal σ µν = diag(σ 11 , σ 22 , σ 33 ) . Pursuant to our background geometry of study (13), it is comfortably deduced that σ 11 = −(σ 22 + σ 33 ) 3 . The propagations of matter parts ( pf / m and em ) and each element of the shear tensor, are then given by 4 Using equation (22) and its first integral, it can easily be indicated that This well-known relation is used frequently in this paper.

Reconstruction of f (R)-gravity
In flat FRW background (i.e. ds 2 = −dt 2 + a 2 (t) FRW (dx 2 + dy 2 + dz 2 ) ), writing the Taylor expansion of the scale factor down as some dimensionless parameters, namely the deceleration, q , jerk, j , snap, s, and others, are appeared and they are defined as [40,41] In (27), a FRW0 and H 0 are the present values of the scale factor and Hubble parameter, respectively. These parameters are a focus of interest because their amounts give important knowledge about the universe. In our anisotropic background geometry, LRS B-I, the relation (27) can be regarded as the Taylor expansion of the average scale factorā = (ab 2 ) 1/3 . Consequently, we may rewrite the dimensionless parameters (28) using the extrinsic curvature Θ . This make our work easy. Furthermore, in order to reconstruct f (R) -gravity models by (22), we 3 These types of shears may usually be defined as where H is the Hubble parameter. 4 Although we present the equations in general form, in examples excluding the last one, the special form of perfect fluid namely pressureless dust matter (i.e. w = 0 ) would be the case of study.
need Θ , not a(t) or b(t) . Hence, we deal with the rate of volume expansion parameter Θ , not the scale factors.
In a spatially homogeneous model the ratio of shear scalar σ to the to expansion scalar Θ is constant: (σ/Θ) = const. . This may compel the following conditions among the directional Hubble parameters H a and H b (H a along x direction while H b along y and z directions) and the expansion scalar Θ : where 1 and 2 are constant, yielding as we expect (because (σ/Θ) = ( ). This physical relation makes our work easy in the calculations. It is important to mention that the special case σ = 0 implies 1 = 2 = (1/3) meaning an isotropic background -FRW. Note that according to (30), the rate of σ and Θ would be equal since (σ/σ) = (Θ/Θ) . Pursuant to observational data, it has been demonstrated in [37,42] that

A constant deceleration parameter (An accelerating universe)
As the first example, let us focus on the acceleration epoch of the universe. This feature of the universe can be determined by the deceleration parameter q -for an accelerating universe q < 0 and for a decelerating universe q > 0 . For the aforementioned purpose, a constant deceleration parameter as is our start point in this part. Writing the deceleration parameter in terms of the extrinsic curvature as and combining it with (31), yields where t 0 is an integration constant which we must take it as an initial time (i.e. at all times of interest t 0 < t) in order to keep expanding phase of the universe (i.e. if t 0 < t then we have Θ > 0 ). Using the average scale factor a , the average Hubble parameterH , and the B-function [27], one has for this case. According to (31) we get According to ref. [27], this era is an accelerated era with non-phantom-like regime property. The limited amounts −1 and 0 correspond to the inflection point namely shifting from a decelerated to an accelerated expansion (i.e. an expansion with constant rate) and de Sitter era/expansion, respectively. It is worthwhile to mention that the left bound, −1 , is equivalent to m = 0 and w eff. = −1/3 and the right bound, 0 , is equivalent to |m| = 1 and w eff. = −1 . Note that we have used the well-known definition in which ρ eff. and p eff. are given by whence we get It is worth mentioning that this relation is also deducible from the relation and equation (33). Using (33), (30), and (26) in (22) one arrives at: where in which The total density in (41), ρ = ρ pf + ρ em , is given by where in which C 1 and C 2 are integration constants. Equations (41) and (45) give the form of f (R) as follows: where C 3 and C 4 are integration constants and The terms l 8 R v and l 9 R n are related to m − part and em − part , respectively. In Fig. 1, the possible powers of R with their generators have been demonstrated for FRW-case. According to (48)-(50), one always has n > v . As is observed, the orders between (1 − √ 3)/2 ≈ −0.366 and (1 + √ 3)/2 ≈ +1.366 are out of access. The inverse powers of R can be generated only by k -parameter. The positive orders of R (here, greater than 1.366)so-called higher orders -can be produced by l -parameter. The matter part can only produce the powers greater than 3/2 while the electromagnetic part can generate greater than +2 . Therefore, practically, the contribution of both m − part and em − part lie on higher orders which can also be generated by l -parameter. It means that the feasible orders of R for a vacuum case and a universe filled with matter/electromagnetic/both matter and electromagnetic, are the same. Note that these discussions are valid only for an ever-accelerating universe. It is important to mention that l , k , and n are affected by anisotropic term, while power v , which comes form matter part, is not affected by it. For anisotropic one, the dashed lines in Fig. 1 shift a little (of order 10 −10 ). Beside all these discussions, it is worthwhile mentioning that there is no way to reach Einstein's theory (i.e. f (R) = R ), and it may back to this point that Einstein's theory does not lead to an ever-accelerating universe. As is clear, from the mathematical viewpoint, all the integration constants and consequently the parameters C 3 , C 4 , l 8 , and l 9 can take any complex value in general. Therefore, let us set them to one 5 ( C 3 = C 4 = l 8 = l 9 = 1 ) for simplicity. It means that the participation amplitude of each term of the obtained form of f (R) in (51) has been normalized to +1 . In order to compare our plots with the results of ref. [1], let us take |m| = 0.5 . As is observed from Fig. 2, the manner of our plots are different than the ones studied in the aforementioned reference either when we consider limited case namely without em − part (see figures 2 and 3 in ref. [1]). The plots are presented for FRW case and note that for the anisotropic one, the plots are shifted so little such that they are not visible (i.e. the total behavior will be unchanged). As we know, from the viability analysis viewpoint, two conditions must be adopted for a f (R) -model: f (R) > 0 (for having a positive effective constant of gravitation) and f (R) > 0 (for the stability of the model). At low curvature, there will be negative and anomalous parts for f and f due to the term R k as k is always negative. Without loss of generality, one may set C 4 = 0 and remove this term, then both validity conditions are satisfied. If one keeps C 4 non-zero, then there will be a lower bound for the evolution range of curvature. For example, in plotting Fig. 2, we kept C 4 , hence there was a lower bound for curvature for the satisfaction of f > 0 as R > 0.2411751621 . For this reason, the related plots have been presented for R ≥ 1 . Keeping C 4 = 0 causes that the growing speed of f (R) be faster than f (R) and f (R) as R increases; for example, at R = 500 the values of f , f , and f are of orders 10 13 , 10 11 , and 10 8 , respectively.

A constant, a variable, and slowly variable jerk parameters
As is observed from (28), another interesting dimensionless parameter for focusing is the jerk parameter. It is 5 Note that for this purpose, C 2 must be taken as a pure imaginary number for having C 2 2 < 0 ; see equations (47)   not hard to show that this parameter in terms of the extrinsic curvature can be written as Three types of jerk values are of physical interest: 1. A constant jerk parameter, j = 1 , which mimics the ΛCDM model; 2. A variable jerk parameter, j = Q(t) , such that the jerk parameter is proportional to the Hubble parameter by an inverse square relation. Hence, Q(t) may be taken as j = Q(t) = λ 2 /H 2 where λ is an arbitrary non-zero real constant and H is the mean Hubble parameter.

A slowly varying jerk parameter.
Since the formulation of two first cases of interest is the same in some parts, hence we first consider these two. The slowly varying case will be considered at the end of this section separately. The solutions to equation (56) for j = 1 and j = Q(t) cases of interest are as follows: The special values n = 2 and n = 3 give j = 1 and j = λ 2 /H 2 , respectively. Note that both (57) and (58) are a solution to each case with the aforementioned conditions. Therefore, for the first case of interest namely j = 1 , the parameter λ, in the above solutions, is a free constant parameter while for the second case namely j = Q(t) it is exactly the constant parameter of relation j = λ 2 /H 2 . For (57) and (58) the effective EoS reads respectively. Therefore, for the constant jerk case ( j = 1 or equivalently n = 2 ), the boundary values will be and for the variable jerk parameter ( j = λ 2 /H 2 or equivalently n = 3 ), we have As is clear from (61) and (63), the solution (57) leads to an unacceptable model because it provides an everaccelerated universe with wrong behavior of w eff. s because both w eff. s stay in negative region and decay from a high value to −1 and consequently it lacks radiation and matter-dominated eras in the past for both jerk types. However, the amounts of w eff. s in (61) and (63), indicate phantom-like regime, but a physical EoS must decrease from positive values (crossing from +1/3 (radiation-dominated era) and 0 (matter-dominated era)) to negative values. Besides these problems, in what follows, it is argued that this type of solution yields an imaginary form for f (R) -gravity which is non-physical. According to (62) and (64) the solution (58) may be called physical as both behaviors of w eff. s are in accordance with a part of the (accelerating) evolution of the observed universe with a non-phantom-like regime. A difference between the two is that w eff.2 = 0 in (62) indicates matter-dominated era that pursuant to w eff.2 = −1 in (62), it is followed by an accelerating expansion, while for the next case of study the start point of EoS is −1/3 which only refers to an acceleration mode of expansion. This difference may be interpreted as a privilege to ΛCDM model in comparison with a decaying jerk model. Note that whether the obtained forms of Θ lead to a physical behavior of EoS or not, it is better we study all these cases because they may yield a form of f (R) which may help for giving a model and solving some problems in other subjects of physics (e.g. Inflation). In other words, if an obtained form of f (R) does not work here because of its nonphysical outcomes, it may be examined in other models which besides f (R) term there is, for example, a scalar field lagrangian and then it may solve some problems. For both sets of solutions, (57)-(58), the Raychaudhuri equation turns out to be where The density in (65) is given by in which where ρ m0 and ρ em0 are constants of integration and For solution (57) ⇔ = +1; For solution (58) ⇔ = −1. In this part, we proceed with a constant jerk parameter. Limpidly, four options are of interest: • ρ m = 0 and ρ em = 0 ; • ρ m = 0 and ρ em = 0 ; • ρ m = 0 and ρ em = 0 ; • ρ m = 0 and ρ em = 0 .
If we keep both densities non-zero, then the analytical solution of equation (65) will be a very complicated case in terms of hypergeometric function. More precisely, besides some hypergeometric functions, there are some complicated analytically unsolvable integrals in terms of hypergeometric functions. The same situation is for the pure electromagnetic case (i.e. ρ m = 0 and ρ em = 0 ). It means that the electromagnetic part leads to some analytically unsolvable integrals. Nonetheless, I think that we may do a thing: "Solving numerically and then fitting the obtained curve with a suitable function in each interval of interest". This, however, provides an approximation function to f (R) but it helps for observing the manner of the evolution of f (R) and giving a model for f (R) -gravity. Between two first options, let us proceed with the generalist case namely ρ m = 0 and ρ em = 0 . It is needless to consider the second option as well since the solving process for both is the same. Solving equation (65) with the conditions numerically via the Runge-Kutta-Fehlberg 4th order method provides three figures 3-5. The figures are presented in the curvature interval 7 through 600 . As is observed, both conditions f (R) > 0 and f (R) > 0 are satisfied in this large interval. According to figures, it seems that polynomial function or some forms of exponential function are good candidates for fitting the curve. On the other hand, according to the plot of f (R) > 0 , the form of f (R) is at least forth order as where c i s are constants. Now, if we fit the form of f (R) in the curvature interval 10 to 100 polynomially with step 1 , we arrive at the form:   Note that the precision of this work depends upon the curvature interval length. Obviously, fitting in a small interval with small step gives a good approximation for the function. We started the interval from 10 instead of 7 because of ignoring some departures observed in the plot of f (R) (see Fig. 5). Our starting point to numerical analysis was R = 7 because (68) and (72) lead to hence, for having physical behavior, the minimum value of curvature is four. According to (74), our approximation type for f (R) looks good for the rest of the interval. It means that at high precision, the form of f (R) would also be a polynomial function.
Two last options namely {ρ m = 0 & ρ em = 0} and {ρ m = 0 & ρ em = 0}, for FRW case have sufficiently been studied in ref. [1]. The solution for anisotropic one is not so different than FRW's solution. Only the constants will vary a little. Hence we do not discuss this case.

The variable jerk case
In this part, we study a variable jerk parameter. Again, like the constant jerk case, here there are four options of interest.
• Pure electromagnetic case ( ρ em = 0 and ρ m = 0): For this case, the solution to the Raychaudhuri equation is found as where A 6 is a constant of integration. This solution holds for both types of extrinsic curvatures in (57)-(58). Note that this solution is obtained for the special case 1 = 2 = 1/3 namely FRW. For the anisotropic case, the basic equation yields an analytically unsolvable integral. Hence, it should be solved numerically. But since the order of anisotropy of the universe is so little, hence its solution will not so different than this solution. Indeed, the solution (76) can also be regarded as a curve fitted function to anisotropic one as well. This situation is in three cases in what follows as well. According to the solution (76), both validity conditions f (R) > 0 and f (R) > 0 are satisfied only by taking A 6 > 0 . The behaviors of f , f and f in the curvature interval [13,35] are presented in Figs. 6-8 with blue color.
• Pure matter case ( ρ em = 0 and ρ m = 0): For this case, the solution to the Raychaudhuri equation for = −1 , which corresponds to (58), is obtained as where A 7 is a constant of integration and e is the base of the natural logarithm, e = 2.718281828 · · · , and " erf " is the error function 6 . For the case = +1 , which corresponds to (57), the corresponded solution yields imaginary form for f (R) which is non-physical, hence we put it aside. As mentioned earlier, the case (57) has further problems as well.
Because the value of error function in every point is in −1 ≤ erf(x) ≤ +1 , hence the last term in (77) is not so strange thing as it can be removed by A 7 . Furthermore, the last term tends to zero as R increases 7 and it means that the effect of this term on the evolution of f (R) is so little such that it may reasonably be ignored. Indeed, the fluctuations and surplusage produced by the last term are practically petty and insignificant. As is clear from (77), by taking R > 12λ 2 , both validity conditions ( f > 0 & f > 0 ) are satisfied. For example, three figures are presented for this case in which the constants have been taken as A 7 = ρ m0 = 1 ; see green plots in Figs. 6-8.
• Both electromagnetic and matter cases (ρ em = 0 and ρ m = 0): 6 The error function is defined for all complex u by erf(u) = 2 √ π u 0 exp(−t 2 )dt . The error function is a smooth function which has a simple zero at u = 0 . 7  In general, this case leads to an integral which is analytically unsolvable. But, by choosing λ = +1 and = −1 , one can arrive at: where A 8 is a constant of integration. This solution is like the sum of the solutions of two cases which obtained above. Limpidly, the error function comes from the matter part not electromagnetic. Also, the root of existing the powers 1/2 and 3/2 backs to matter part while the electromagnetic part produces +2 instead. However, the behavior of (78) depends upon the selection of constant parameters, but obviously, in general, both validity conditions for (78) are satisfied by taking R > 12 . Taking start point R = 13 and setting A 8 = 1 and ρ m0 = ρ em0 = 1/2 , three plots have been presented in the interval 13 through 35 ; see red plots in Figs. 6-8.
where A 9 is a constant of integration. This solution is arguable from the last three solutions as the term which does not depend upon the densities is this term. Hence, it is the effect of the vacuum case. Clearly, both validity conditions are satisfied only by taking A 9 > 0 . Unlike three last cases, here, there is no lowers bound for starting the physical values of curvature according to (79).
The related plots to this case have been demonstrated by orange color in Figs. 6-8. According to all solutions obtained for j = Q(t) and their Taylor expansions 8 , all four f (R) s of this case obviously reduce to Einstein's theory at so low curvature. Furthermore, at low curvature, f , f , and f of four options in view of their amounts, satisfy this relation: em 9 > m + em 10 > vacuum 11 > m 12 (see Figs. [6][7][8]. At high curvature, the arrangement will be changed as em > vacuum > m + em > m . It means that adding electromagnetic part causes that the amounts of f , f , and f increase, while adding typical matter/perfect fluid implies that their amounts decrease than vacuum case. In all curvature of interest, f , f , and f of pure electromagnetic and pure matter cases have the highest and lowest values, respectively. At high curvature, f , f , and f of vacuum case are greater than the case which contains both matter and electromagnetic while at low curvature, the layout is changed.

A slowly varying jerk parameter
In this part, a scenario where the jerk parameter is a slowly varying function of redshift z , viz, where η 1 is a small constant parameter and F (z) is a slowly varying function of the redshift, is considered. Since j(z) varies slowly with respect to z , hence it is a good approximate that we take F (z) ≈ F 0 + F 1 z ( F 0 and F 1 are constant). Therefore, one may easily getΘ 8 We know that erf(u) = 2 k!(2k+1) . 9 Pure electromagnetic case 10 Both electromagnetic and matter cases 11 Vacuum case 12 Pure matter case where η 3 and η 4 are integration constants and Obviously, the powers related to vacuum part namely ξ 1 and ξ 2 can take any constant, but it is better we take them positive for the satisfaction of validity conditions. The participation of electromagnetic and perfect fluid are as R 4|β1| and R 2(1+w) , respectively. Again like previous cases of study, we observe that unlike the matter part term, the power of curvature arisen from the electromagnetic part is affected by the anisotropic background. Note that |β 1 | ≈ 2/3 hence the electromagnetic part for this case produces R 8/3 ( 8/3 is the exact value of FRW background). The special cases namely pressureless dust matter (w = 0 ) generates R 2 which is less than the electromagnetic case. Therefore, pursuant to the previous examples and this one, it seems that the electromagnetic case in most cases of interest, generates the higher powers of curvature than perfect fluid/matter case. It is interesting to note that for the dark energy ( w = −1 ) the power 2(1 + w) would be zero, hence it does not produce any power for the curvature. Both perfect fluid and electromagnetic parts satisfy the validity conditions separately. And as a final point, we mention that Einstein's general relativity theory can be recovered by the terms of the vacuum part or by a matter with the EoS w = (−1/2) .

On the enhancement of the method
The presented method has this potential to be applied further. The first generalization is that the inverse road of the method may be adopted to examine a given theory of f (R) according to its emerged outcomes (i.e. behaviors of shear and Hubble, etc.) and comparing them with observational data. But, it seems that in most cases of interest, this application is numerically feasible, not analytically.
According to observational data, the treatment of Hubble, EoS, and etc. are clear. On the other hand, we have Raychaudhuri equation which gives us the form of f (R) . Hence, by applying a curve-fitting method to obtained curve from numerical methods, it is feasible to arrive at some forms for f (R) at different stages of the evolution of the universe. This is the second generalization.
Because these tasks are beyond the scope of this paper, hence I do not give an example but pursuant to the first suggestion, I tried for some given forms of f (R) and found that they are only numerically doable.

Conclusions
Utilizing Raychaudhuri-based reconstruction strategy suggested by Choudhury et al. [1], the reconstruction of anisotropic Einstein-Maxwell equation in 1 + 3 covariant formalism of f (R) -gravity was investigated. The matter part of the problem was assumed to be a non-interacting combination of a perfect fluid and an electromagnetic field. The model has been reconstructed in four interesting modes of evolution. In summary, some of our findings to these modes were as follows: 1. A constant deceleration parameter (An accelerating universe): The obtained form for f (R) was as: The range of powers was given in Fig. 1. It has been concluded that Einstein's theory does not emerge from this form because of the domains of powers and it is due to the fact that his theory does not give an everaccelerating universe. The terms l 8 R v and l 9 R n in the above form come from matter and electromagnetic parts, respectively. Under the conditions of the problem, one always has n > v and both contribute at the higher orders of curvature which are also reachable via vacuum case.
2. The constant jerk case j = 1 mimicking ΛCDM model: Generally, this case is analytically unsolvable. Hence, we proceed using the Runge-Kutta-Fehlberg 4th order method and (functional) curve-fitting method. The outcome for the generalist case up to fifth order polynomial in some curvature interval was as: The obtained forms to f (R) were respectively as follows: Pursuant to these forms, reconstruction via variable jerk led to an exponential function of f (R) , exp(R−R 0 ) , and the contribution of matter and electromagnetic parts appeared as respectively. All these interesting obtained forms at low curvature tend to Einstein's theory of gravity.
The first two terms come from vacuum part while the third and fourth terms arise from the electromagnetic and perfect fluid, respectively. Hence, the participation of the electromagnetic part for FRW is as R 8/3 . For the pressureless dust matter, the power of curvature will be +2 and for the dark energy, the power of curvature is zero. The validity conditions for all f (R) s obtained from four modes of evolution were satisfied (in some cases entirely and in others at special intervals or under specific conditions).
There is a interesting common property among all cases studied in this paper: Unlike the perfect fluid/matter part, the power of curvature produced by the electromagnetic part is affected by anisotropic property of background. Furthermore, the power of curvature supplied by electromagnetic part is higher than matter/perfect fluid part. For example, for FRW case, some of our findings were as follows: • For constant deceleration parameter: P em > 2 and P m > 3/2 ; • For variable jerk parameter with time: P em = 2 and P m = 3/2 & 1/2 ; • For slowly varying jerk parameter with redshift: P em = 8/3 and P m = 2 , where P em and P m refer to the powers of curvatures of the electromagnetic and matter parts, respectively.
Finally, some discussions about the enhancement of the method were done.