Analytic charged BHs in $f(\mathcal{R})$ gravity

In this article, we seek exact charged spherically symmetric black holes (BHs) with considering $f(\mathcal{R})$ gravitational theory. These BHs are characterized by convolution and error functions. Those two functions depend on a constant of integration which is responsible to make such a solution deviate from the Einstein general relativity (GR). The error function which constitutes the charge potential of the Maxwell field depends on the constant of integration and when this constant is vanishing we can not reproduce the Reissner-Nordstr\"om BH in the lower order of $f(\mathcal{R})$. This means that we can not reproduce Reissner-Nordstr\"om BH in lower-order-curvature theory, i.e., in GR limit $f(\mathcal{R})=\mathcal{R}$, we can not get the well known charged BH. We study the physical properties of these BHs and show that it is asymptotically approached as a flat spacetime or approach AdS/dS spacetime. Also, we calculate the invariants of the BHS and show that the singularities are milder than those of BH's of GR. Additionally, we derive the stability condition through the use of geodesic deviation. Moreover, we study the thermodynamics of our BH and investigate the impact of the higher-order-curvature theory. Finally, we show that all the BHs are stable and have radial speed equal to one through the use of odd-type mode.


I. INTRODUCTION
Despite the simplicity, elegantly and powerfully of GR, which is considered as a classical theory, and despite its successes in the prediction of the shadow of BH, which confirmed by horizon telescope [1][2][3], and its success with observation, like the bending of light, perihelion precession of mercury, the gravitational redshift of radiation from distant stars [4][5][6], it has unsolved problems. Among these issues, is the problem of dark energy and dark matter which confirmed by recent observation [7][8][9][10][11], the problem of singularities [12][13][14], the problem of quantum gravity [15,16]. All these issues make physicists demand another theory that is eligible and able to overcome such problems and permit GR in the lower energy. Therefore, a satiate modified gravitational theory has been constructed that can overcome the issues that GR cannot solve. Later, physicists put constraints that any modified gravity must bypass like: its consistency with the solar system, free from ghost modes, able to describe the observations in which GR could not solve, and must not create a fifth force in local physics. There are many modified theories that satisfy such benchmark and are categorized as: a) Theories in which their Lagrangian involve higher-order curvature like f (R), [17,18], Lanczos-Lovelock models [19][20][21][22][23] etc.
In this study, we take f (R) gravitational theory which is considered as the excellent modification of the Einstein GR at the most recent time. Among many things that make f (R) gravity be the best modification of GR is its simplified Lagrangian which contains an arbitrary function whose lowest order coincides with Einstein's theory. The f (R) gravity is an elegant theory to describe an expanded accelerated universe which is supported from recent observations [11,34]. Moreover, it is supposed that a special form of f (R) can mimic some curvature modifications that come from the quantum theory of gravity. The justifications of invented modified theories of GR are to construct complete systems of gravitational models to cure the issues appearing in GR. These systems are considered in various topics of gravity like astrophysics, cosmology and the output of these studies construct as a suitable extension of the Einstein GR [35][36][37]. Recently, the f (R) gravity theory gains much attentions which generalize the gravitational Lagrangian of GR [38][39][40][41]. There are two methods to obtain the equation of motions of f (R) gravitational theory, the first through the use of metric structure and the other through the use of the Palatini method in which the connection and metric are taken as independent variables [42]. There are many viable applications of f (R) in the domain of cosmology as well as astrophysics [43][44][45][46] and in a specific case, f (R) can reduce to GR case [47].
There are many spherically symmetric BH solutions derived in f (R) [48][49][50][51][52]. Moreover, novel spherically symmetric BH solutions have been obtained in the context of f (R) [53][54][55]. It is well known that higher-order curvature terms play an ingredient role in strong gravitational background in local objects of f (R). In this context, many physicists concentrate on the study of spherically symmetric BHs [56][57][58][59][60][61][62][63]. This exposition aims to derive original spherically symmetric charged BH solutions in f (R) gravity devoid of any constraints on the form of f (R) nor the Ricci scalar and study their contingent physics. In the frame of f (R)-gravity, a Lagrangian approach has been developed to study dynamics of spherically symmetric metrics. The Euler-Lagrange equations are obtained and solved in the case of constant curvature R = R 0 recovering the standard Schwarzschild-de Sitter solution of GR [64].
The paper is structured as: In Sec. II, we give a brief construction of f (R) gravity. In Sec. III, we employ the charged equation of motions of f (R) to a spherically symmetric spacetime having two unknown metric potentials and derive a system of non-linear differential equations. These systems have four unknown functions, two of the metric potentials, one from the gauge potential of the charge and the derivative of f (R). To be able to solve this system and to put it in a closed system, we assume a form of the derivative of f (R), which depends on a constant of integration and solve the system. The BH solution is characterized by convolution and error functions and they become constant functions when the constant of integration that appears in the derivative function of f (R) is equal to zero and in that case, we have a GR BHs. Therefore, the convolution and error functions appear like the effect of higher-order curvature that characterizes f (R) gravity. Moreover, we study the physics of these BHs by giving the asymptote form of the convolution and error functions up to certain order and showing the effect of the higher-order curvature. We examine the singularities of the BHs by calculating the invariants of the Kretschmann scalar, the Ricci tensor square, and the Ricci scalar and investigate the effect of higher-order curvature on these invariants by showing that the singularity of our BHs is milder than the GR BHs. Also in Sec. III, we calculate the geodesic deviation and give the condition of stability. In Sec. IV, we show the validity of the first law of thermodynamics and also calculate some thermodynamic quantities like the Hawking temperature, entropy, the Gibbs free energy, and quasi-local energy of these BHs. We reserve the final section for the discussion and conclusion of this study.

II. FUNDAMENTALS OF f (R) GRAVITATIONAL THEORY
We are going to study a 4-dimensional Lagrangian of f (R) with f being an arbitrary function. The Lagrangian of f (R) is given by (cf. [17,40,[65][66][67][68][69][70]): Here κ and g are the Newtonian gravitational constant and the determinant of the metric. The Lagrangian of electromagnetic field L em is defined as L em = F where F = dξ and ξ = ξ β dx β is the electromagnetic Maxwell gauge potential 1-form [71].
Making the variations principle w.r.t. the metric tensor g αβ to the Lagrangian (1) we get the equation of motions of f (R) gravitational theory in the form [72] with ✷ being the d'Alembertian operator, f R = df dR , and T (em) ν µ is the energy-momentum tensor of the Maxwell field defined as where F = F µν F µν . Moreover, the variation of Eq. (1) w.r.t. the gauge potential 1-form, ξ µ , gives: The trace of the field equations (2), takes the form: From Eq. (5) one can obtain f (R) in the form: Using Eq. (6) in Eq. (2) we get [73] We will apply the field equations (5), (7), and (3) to a spherically symmetric spacetime with two unequal unknown functions in the next section.

III. SPHERICALLY SYMMETRIC BLACK HOLE SOLUTIONS
Let us take the following line-element which describes a spherically symmetric a spacetime with µ(r) and ν(r) which are two functions of the radial coordinate r to be determined from the equation of motions. From Eq. (8), the Ricci scalar takes the form: where µ ≡ µ(r), ν ≡ ν(r), µ ′ = dµ dr , µ ′′ = d2µ dr 2 , and ν ′ = dν dr . Applying the equation of motions (4), (5) and (7) to the line-element (8) using Eq. (9), we get the following non-linear differential equations: where α(r) = µ(r) ν(r) and The non-vanishing components of the Maxwell field has the form When we exclude the trace part, Eqs. (10) can have the form: Using Eqs. (12) and (13), i.e., [(12) minus (13)], we get Moreover, Eqs. (12) and (14), i.e., [(12) plus (14)] give: Careful check can show that Eqs. (16) and (15) are coincide and therefore we get two equations from Eqs. (12), (13) and (14) which are independent. Using the above information, we can say that Eq. (12) is equal to minus Eq. (13) minus two time Eq. (14). Thus, Eqs. (12) and Eq. (16) can be chosen as independent equations. Due to the fact that we have four unknown functions ν, α, ξ, and F , we cannot determine one function. Let us discuss some special cases of Eq. (15) or Eq. (16). we can get the Reissner-Nordström solution by supposing, Assuming ν = 0 we can get from Eq. (15) Using Eq. (18) in Eq. (12) we get Assuming F 3 = 0, Eq. (19) gives From Eq. (20) after using Eq. (11) we get the following solution where ν 0 ν 1 and ν 2 are integration constants. Equation (21) is the well-known Reissner-Nordström-(anti-)de Sitter space-time. Moreover, we study the case F 2 = 0 in which Eq. (19) gives, The solution of Eq. (22) together with Eq. (11) have the following solution whereν 0 ,ν 1 andν 2 are constants. Equation (23) corresponds to the solution derived before in [53,55]. For the general case, i.e., when F 2 and F 3 are not vanishing and for small r then F 2 term in Eq. (19) will be dominates and the solution must has the form given in Eq. (21). When r is large then F 3 term in Eq. (19) will be dominated and the solution has the form given by Eq. (23). Thus, there should be a solution that connects the two solutions given by Eqs. (21), for a small r region, and (23) for large r region.
Because we are studying a spherical symmetry case, we assume f (R) = f (r). The above system of differential equations, Eqs. (10) and (11), are five non-linear differential equations in four unknowns, F 1 , ν, α, and ξ. Thus, we cannot solve such a system unless we adjust the unknown with the number of differential equations. The only way to solve the above system is to exclude the trace equation, i.e., we must not take into account the differential equation I from Eq. (10) and in that case, the solution of the system becomes: where a 0 , a 1 , a 2 , and a 3 are constants. Equation (24) shows that when a 1 = 0, we return to GR case. Here, Finally the erf function is the error function which is defined by Using Eq. (24) in the trace equation, i.e., the fourth equation of Eq. (10) we get f (r) in the form 1 The HeunC function is the solution of the Heun Confluent equation which is defined as The solution of the above differential equation defined HeunC(α, β, γ, δ, η, r) for more details, interested readers can check. The He-unCPrime is the derivative of the Heun Confluent function.
Using Eq. (26) in Eq. (9), we get Equations (24), (26), and (27) show that when a 1 = 0, we get Equation (28) indicates that when F 1 (r) = 1, we get f (R) = R which gives µ(r) = ν(r) = 1 + a2 r + a 3 r 2 . All the above informations guarantees that when a 1 = 0 we get the GR black holes which makes the constant a 3 acts as a cosmological constant 2 . Equation (28) ensures that in higher-order curvature theory we cannot reproduce the ordinary Reissner-Nordström. This is because that the charged solution in this class of modified theory is related to the higher order curvature, i.e., depends on a 1 and if a 1 = 0, we will not get the Reissner-Nordström solution of GR.
A. Physical properties of the black hole (24) In this section, we try to understand the physical properties of the black hole (24). The asymptote behaviors of the metric potentials, µ(r) and ν(r) in Eq. (24), by assuming a 3 = 0, have the form In Eq. (29), we have put a 2 = −2M . Using Eq. (29) in Eq. (8), we get The line element (30) is asymptotically approaching a flat spacetime and does not coincide with the Reissner-Nordström spacetime due to the contribution of the extra terms that come mainly from the constant parameter a 1 whose source is the higher-order curvature terms of f (R) gravitational theory. Equation (30) shows in a clear way that in the higher-order curvature gravity, one can get a spacetime different from the Reissner-Nordström and when the constant a 1 = 0, we can recover the Reissner-Nordström metric. As one can check easily that when these extra terms equal zero one can smoothly return to the Schwarzschild spacetime. In conclusion, we can say that in higher-order curvature gravity, we can get a charged spacetime that is different from the Reissner-Nordström and cannot reduce to the Reissner-Nordström in the lower order of f (R) = R. Now we going to assume that the constant a 3 = 0 and get the metric potentials in the form 2 Note that when a 1 = 0, we get H = H 1 = HeunC 3 2 , 3 2 , 0, 3 8 , 9 8 , 0 = HeunC 3 2 , − 3 2 , 0, 3 8 , 9 8 , 0 = 1, H 2 = HeunCPrime 3 2 , 3 2 , 0, 3 8 , 9 8 , − a 1 r 2 = 0 and H 3 = HeunCPrime 3 2 , − 3 2 , 0, 3 8 , 9 8 , − a 1 where Λ eff = a 3 . Uses of Eq. (31) in (8) gives The line element (30) is asymptotically approaches the AdS/dS spacetime according to the sign of Λ eff which depends mainly on the constant a 3 . Now we are going to use Eq. (32) in Eq. (9) and get We have neglected the other two roots of Eq. (33) because they are imaginary. From Eq. (33), it is clear that when the constant a 3 = 0 we get a non-vanishing value of the Ricci scalar because of the higher-order curvature. When a 1 = 0 we get a vanishing value of the Ricci scalar which corresponds to GR black hole. The asymptote form of the function f (r) given by Eq. (26) becomes Using second equation of (33) in (34), we get where C i , i = 1, · · · , 4 are lengthy constants whose values depend on a 1 and a 3 . Using Eq. (31), we get the invariants of solution (24) as 3 Here (R µνρσ R µνρσ , R µν R µν , R) represent the Kretschmann scalar, the Ricci tensor square, the Ricci scalar, respectively and all of them have a true singularity at r = 0. It is important to note that the constant a 1 is the source of the deviation of the above results of charged solution from the Reissner-Nordström spacetime of GR whose invariants have the form (R µνρσ R µνρσ , R µν R µν , R) = 24Λ 2 + 48M 2 r 6 , 36Λ 2 + 4q 4 r 8 , 12Λ . Equation (36) indicates that the leading term of the invariants (R µνρσ R µνρσ , R µν R µν , R) is 1 r 2 , 1 r 2 , 1 r 2 which is different from the Reissner-Nordström black hole whose leading terms of the Kretschmann and the Ricci tensor squared are 1 r 6 , 1 r 8 . Therefore, Eq. (36) indicates that the singularity of the Kretschmann scalar and the Ricci tensor squared are much softer than those of Reissner-Nordström black hole of GR.

B. Stability of the black holes using geodesic deviation
The geodesic equations of a test particle in the gravitational field are given by where s represents the affine connection parameter. The geodesic deviation equations have the form [74] where ǫ ρ is the 4-vector deviation. Introducing (37) and (38) into (32), we get and for the geodesic deviation, the black hole spacetime (32) gives where µ(r) and ν(r) are defined in Eq. (32) and ′ means derivative w.r.t. the radial coordinate r. Using the circular orbit we get Equations (40) can be rewritten as The second equation of (43) corresponds to a simple harmonic motion, which means that there is stability on the plane θ = π/2. Assuming the remaining equations of (43) to have solutions in the form where ζ 1 , ζ 2 , and ζ 3 are constants and ϕ is an unknown variable. Using Eq. (44) into (43), one can get the stability condition for a static spherically symmetric charged black hole in the form Equation (45) can has the following solution We depicted Eq. (46) in Fig. 1 for particular values of the model. The case Λ eff = a 3 = 0 is drawn in Fig. 1 (a) and the case Λ eff = a 3 = 0 is drawn in Fig. 1 (b). These two figures exhibit the regions where the black holes are stable and not stable.  (29) and (31).
The black hole (29) is characterized by the mass M and the parameter a 1 and when the parameter a 1 is vanishing we get the Schwarzschild black hole 5 which corresponds to GR. To find the horizons of the black hole, (29), we put µ(r) = 0 and get four lengthy real positive roots. The metric potentials of the black hole (29) are drawn in Fig. 2 (a). From Fig. 2 (a), we can easy see the two positive horizons of the metric potentials µ(r). It is easy to check that the degenerate horizon for the metric potential µ(r) is happened for a specific values of (a 1 , M, r) ≡ (−0.043, 0.5, 0.5), respectively, which correspond to the Nariai black hole. The degenerate behavior is shown in Fig. 2 (b). Fig. 2 (b) shows that the horizon r 1 increases as a 1 increases while r 2 decreases as a 1 decreases. As we observe from Fig. 2 (b) that as a 1 increases and M decrease i,e, a 1 > M , we enter in a parameter region where there is no horizon, and thus the central singularity is a naked singularity. Using Eq. (47), we get the Hawking temperature in the form  (29) ; (c) typical behaviour of the horizon temperature, (51), which shows that T1 has a positive decreasing value while T2 has a negative increasing value; (d) typical behaviour of the horizon entropy, (52), which shows that ψ (1,2) always has positive values, ψ1 increases with M wile ψ2 has positive decreasing value; (e) typical behaviour of the quasi-local energy, (53), that indicates that E (1,2) has positive increasing values and also E1 > E2; (f) typical behaviour of the horizon Gibbs free energy which shows that G1 is positive while G2 starts with positive value and as M increase it becomes negative.
The behavior of the Hawking temperature given by Eq. (51) is drawn in Fig. 2 (c) which shows that T 1 > T 2 . As Fig. 2 (c) shows that the T 1 has a decreasing positive temperature while T 2 has an increasing negative temperature.
Using Eq. (48) we get the entropy of black hole (29) in the form The behavior of the entropy is given in Fig. 1 (d) which shows an increasing value for ψ 1 and decreasing value for ψ 2 . Using Eq. (49) we calculate the quasi-local energy and get The behavior of the quasi-local energies are shown in Fig. 1 (e) which also shows positive increasing values for E (1,2) and also shows E 1 > E 2 . Finally, we use Eqs. (51), (52) and (53) in Eq. (50) to calculate the Gibbs free energies. The behavior of these free energies are shown in Fig. 1 (f) which shows positive increasing value for G 1 and G 2 start with positive decreasing value then it becomes negative value.  A. Thermodynamics of the black hole (31) that has asymptote flat AdS/dS In this subsection, we are going to study the black hole (31) which is characterized by the mass of the black hole M , the parameter a 1 and a positive cosmological effective constant 6 . The metric potentials of the black hole (54) are drawn in Fig. 3 (a). From Fig. 3 (a), we can easily see the two horizons of the metric potentials µ(r) and ν(r). To find the horizons of this black hole, (54), we put µ(r) = 0 in Eq. (54). This gives four roots two of them are real and the others are imaginary. These real roots are lengthy however, their behaviors are drawn in Fig. 3 (b). It is easy to check that the degenerate horizon for the metric potential µ(r) given by Eq. (54) is happened for a specific values for (a 1 , M, r, Λ eff ) ≡ (−0.035, 0.5, 0.3881986644, 1), respectively which corresponds to the Nariai black hole. The degenerate behavior is shown is Fig. 3 (b) which shows that the horizon r 1 increases with a 1 while r 2 decreases. Using Eq. (47) we draw the behavior of the Hawking temperatures in Fig. 3 (c) which shows that T 1 > T 2 . As Fig. 3 (c) shows that the T 1 has decreasing positive temperature while T 2 has increasing negative temperature.
Using Eq. (48), we draw the entropy in Fig. 3 (d) which shows increasing value for ψ 1 and decreasing value for ψ 2 . Using Eq. (49) we calculate the quasi-local energy that is drawn in Fig. 3 (e) which also shows positive increasing value for E (1,2) . Finally, we use Eqs. (51), (52) and (53) in Eq. (50) to calculate the Gibbs free energies. The behavior of these free energies are shown in Fig. 3 (f) which shows negative increasing value for G 1 while G 2 starts with negative 6 The metric potential when Λ eff. = a 3 take the form value and as M increases it becomes positive. By the same procedure we can study the case dS, i.e., Λ eff a 3 < 0 and show the behaviour of their thermodynamical physical quantities.
B. First law of thermodynamics of the BHs (29) and (31) It is important for any black hole to check its validity of the first law of thermodynamics. Therefore, we will use the formula of this law in the frame of f (R) as [82] with E, ψ, T , P , and V are the quasi-local energy, the Hawking entropy, the Hawking temperature and the radial component of the stress-energy tensor that serves as a thermodynamic pressure, P = T r r | ± , and the geometric volume respectively. The pressure, in the context of f (R) gravitational theory, is figured as [82] For the flat spacetime (29) if we neglect O 1 r 3 and O 1 r 4 to make the calculations more applicable, we get 7 Using (58) in (56), we can prove that the first law of the flat spacetime is verified. Repeating the same procedure for the AdS/dS, we get Substituting the form of r into the thermodynamical quantities we get If we use Eq. (59) in (56), we can verify the first law of thermodynamics for the BH (29). If one repeats the same procedure for the BH (31), one can verify the first law of thermodynamics provided that we neglect all the quantities containing b to make the calculation easier to carry out. 7 When we neglect the terms of order O 1 and when A(r) = 0 we get two roots.

V. DISCUSSION AND CONCLUSIONS
It is well-known that spherically symmetric spacetime is an essential spacetime in GR for many reasons: First because it is the easiest geometry that one can start to use in the study of GR. In this exposition, we tackle the problem of applying the charged field equations of f (R) to a spherically symmetric spacetime with two unknown functions. We derive a highly nonlinear five differential equations in four unknowns, two of the metric potential, the third is the derivative of f (R), and the fourth is the gauge potential of the Maxwell field. We study some limits of this system of differential equations to test its compatibility with the previous studies and got consistent results of BHs derived for constant and non-constant Ricci scalar. Then we turn our attention to the general case and solve this system by assuming a specific form of the derivative of f (R). This assumption involves one constant that when it equals zero we got the GR theory limit. We derive the form of the metric potential and the gauge potential of the Maxwell field and show they depend on the convolution and error functions respectively. Under some constraints, if these functions are set zero, we got to the Schwarzschild BH of GR which means that in the limit of higher-order curvature theory, we cannot generate the Reissner-Nordström BH. This gives us an indication that the effect of the higher-order curvature terms restored in the convolution and error functions.
To gain some physics of these BHs we studied their asymptotes and show that we have a constant of integration that plays an important role in the asymptote, i.e., when this constant is vanishing/non-vanishing, we got two different asymptotes of the BHs, a flat and AdS/dS spacetimes respectively. Moreover, we studied the invariants of these BHs and showed that their singularities are milder than those of GR BH. By the use of the geodesic deviation, we derived the conditions of stability of those BHs and showed their stable areas as shown in Figs. 1 (a) and (b). Also, we have calculated some thermodynamics quantities like the Hawking temperature, entropy, quasi-local energy, and the Gibbs free energy to know the behavior of those BHs analytic and graphically and show in detail that such BHs satisfy the first law of thermodynamics. Finally, if we use the results of odd perturbation presented in [53], we can show easily that the BHs derived in this study are stable.
To conclude, we have derived new charged BHs in the context of f (R) that their Ricci scalars are not constant. The originality of those BHs comes from the constant that involves the assumption of the derivative function of f (R). If we follow the same procedure but with another assumption of the derivative function of f (R), we can get new BHs but this needs more study because of the complication of the system. This will be done elsewhere.