Naked singularities in quadratic f(R) gravity
Enrique Hernandéz-Lorenzo Christian F. Steinwachs
gr-qc hep-th




FR-PHENO-2020-004



christian.steinwachs@physik.uni-freiburg.de

Physikalisches Institut, Albert-Ludwigs-Universität Freiburg,

Hermann-Herder-Straße 3, 79104 Freiburg, Germany





We find new static, spherically symmetric, and asymptotically flat vacuum solutions without horizon in Starobinsky's quadratic f(R) gravity.
We systematically classify these solutions by an asymptotic analysis around the origin and find seven different integer Frobenius families.
We numerically solve the exact equations of motion by a double-shooting method and specify boundary conditions by matching the numerical solution to the analytic solution of the linearised field equations in the weak field regime.
We find that all integer Frobenius families can be connected to asymptotically flat solutions and trace out lines in the parameter space, allowing to ultimately relate all free parameters to the total mass at infinity.   





Naked singularities in quadratic f(R) gravity
    Christian F. Steinwachs
    Received November xx, 2019; accepted xxxx xx, 2019
======================================================






§ INTRODUCTION



Theories incorporating geometric modifications of General Relativity (GR) are highly relevant in the context of quantum gravity and cosmology.
In particular, f(R) modifications of GR give rise to phenomenologically vital cosmological models for the early and late time acceleration of the Universe <cit.>, of which Starobinsky's quadratic  f(R)=R+R^2 model of inflation <cit.> is the most relevant one and favoured by CMB data <cit.>.
The higher derivative structure of f(R) gravity leads to an additional scalar degree of freedom, the scalaron, made manifest in the classically equivalent scalar-tensor representation.
In <cit.>, it was shown that this classical equivalence also extends to the quantum level in a similar way as the on-shell quantum equivalence between different field parametrizations in scalar-tensor theories found in <cit.>.
In contrast to higher-derivative theories such as Quadratic Gravity (QG), the degeneracy structure of f(R) gravity ensures the absence of the Ostrogradski instability <cit.>. 

In order to investigate the construction of modified theories of gravity, it is important to compare the structure of solutions to that of GR. In general, as for GR, analytical solutions are only available for a high degree of symmetry.
However, in contrast to GR, in general no uniqueness theorems are available in higher derivative theories of gravity, see e.g. <cit.> and references therein.
Nevertheless, in quadratic gravity (QG) there is a “trace no-hair theorem” <cit.>, which implies that static asymptotically flat vacuum solutions with a horizon must have vanishing Ricci scalar R.
Therefore, in the search for black hole solutions, the R^2 term does not contribute to the equations of motion and the black hole analysis in QG effectively reduces to that in Einstein-Weyl gravity (EWG) <cit.>.
In contrast, we do not search for black hole solutions, but instead focus on solutions without a horizon.
Solutions without a horizon but a singular geometry feature naked singularities and might have interesting implications for the quantum information paradox <cit.> and might serve as candidates for dark matter <cit.>.

In this letter, we systematically investigate and classify static, spherically symmetric and asymptotically flat vacuum solutions in Starobinsky's quadratic f(R) theory by a combined analytic and numerical approach.
We first analytically perform the asymptotic analyses in the near origin regime and in the linearised weak field regime and then connect the two asymptotic regimes by numerically solving the full non-linear field equations by a “double-shooting” algorithm.
In particular, we do not make any additional simplifying assumptions about the curvature or the metric functions in the spherically symmetric ansatz of the line element, which allows us to screen the full space of solutions.





§ STAROBINSKY'S QUADRATIC F(R) MODEL




The f(R) action in four spacetime dimensions reads

    S[g]=∫d^4x √(-g) f(R).

Derivatives of f with respect to R are denoted by f_n:=∂^nf(R)/∂ R^n and the equations of motion (EOM) of f(R) gravity read ℰ^μν=0 with the extremal tensor

    ℰ^μν:=   (g^μν-∇^μ∇^ν) f_1+R^μνf_1-1/2g^μνf.

The trace of the extremal tensor is defined as

    ℰ:= g_μνℰ^μν=3 f_1+Rf_1-2f.

The invariance of the action (<ref>) under diffeomorphisms implies the generalized Bianchi identity

    ∇_μℰ^μν=0.

The fourth-order nature of the field equations implies that in addition to the massless spin-two graviton, f(R) gravity propagates the massive spin-zero scalaron.
We consider the Starobinsky model defined by (<ref>) with 

    f(R)=M_P^2/2(R+R^2/6M_S^2).

The two mass scales (in natural units c=ħ=1) are set by the reduced Planck mass M_P=1/√(8π G_N) and the scalaron mass M_S^2≥0.
The inequality ensures that the scalaron is not a tachyon and M_S=0 corresponds to the induced gravity limit M_P→0 of (<ref>) resulting in the scale invariant R^2 model.
In case of a scalaron-driven inflationary phase of the early universe, M_S is constrained by Cosmic Microwave Background (CMB) data <cit.>,

    M_S≈ 10^-5M_P≈ 10^13GeV.
 




§ STATIC SPHERICALLY SYMMETRIC METRIC ANSATZ




In Schwarzschild coordinates (t,r,ϑ,ϕ), the static and spherically symmetric line element is parametrized as

    ds^2=-B(r)dt^2+A(r)dr^2+r^2dΩ^2_(2).

Here, dΩ^2_(2)=dϑ^2+sin^2ϑdϕ^2 is the line element of the unit two-sphere and A(r) and B(r) are arbitrary functions of the radius r.
Only two of the non-zero components ℰ^tt, ℰ^rr, ℰ^ϑϑ and ℰ^ϕϕ are independent due to the redundancy  ℰ^ϕϕ=sin^2ϑℰ^ϑϑ and the generalized Bianchi identity (<ref>).
We work with the EOMs ℰ^tt=0 and ℰ^rr=0, which have the functional dependencies

    ℰ^tt(B_4,B_3,B_2,B_1,B,A_3, A_2,A_1,A, r;M_S)=   0,
    ℰ^rr(B_3,B_2,B_1,B,A_2,A_1,A,r;M_S)=   0.

The nth derivatives of B and A with respect to r are denoted by B_n and A_n. The explicit expressions for ℰ^tt and ℰ^rr are provided in (<ref>) and (<ref>).
In the form (<ref>), (<ref>), the differential order of this system is not transparent.
From (<ref>), one might expect that the differential order of (<ref>) and (<ref>) is four.
While correct in the present case, the general situation is more complicated, as (<ref>), (<ref>) is not a system of ordinary differential equations (ODE) but a system of differential-algebraic equations (DAE), which involves algebraic constraints.
The DAE structure is made manifest be rewriting  (<ref>), (<ref>) as system of first order equations. 
The differential order D of a DAE with M first-order ODEs, N constraints and differential index I is <cit.>,

    D=M-(I-1)N.

The differential index I is defined as the minimum number of differentiations required to write the  DAEs as ODEs.
For the EOMs  (<ref>), (<ref>), we find M=7, N=1 and I=4.
Hence, D=4 and a general solution of (<ref>), (<ref>) will generically depend on four integration constants.

Following the treatment of <cit.> in the context of QG, for the numerical analysis it is more convenient to work with a system of two second-order EOMs

    E^tt(B_2,B_1,B, A_1, A,r;M_S)=   0,
    
    E^rr(B_2,B_1,B, A_2,A_1,A,r;M_S  )=   0.

The E^tt and E^rr equations are defined as[Modulo possible singularities due to zeros in the denominators of X, Y and Z, the systems (<ref>), (<ref>) and  (<ref>), (<ref>) are equivalent.]   

    E^tt:=   ℰ^tt-Xℰ^rr-Y∂_rℰ^rr,
    
    E^rr:=   ℰ^rr-Z∂_rE^tt.

With a_1:=r∂_rln A, b_1:=r∂_rln B, b_2:=r∂_rb_1, we have

    X:=   -B/A[1+24-2(2+a_1)(4+b_1)/(4+b_1)^2],
    
    Y:=   -2Br[A(4+b_1)]^-1,
    
    Z:=   [A (b_1+4)^3 r]{6 B [2 b_1 (-A M_S^2
    	r^2+b_1+b_2+2)..
       ..+4 (b_2-2 A M_S^2 r^2)-a_1
    	(b_1+4) b_1+b_1^3]}^-1.
 
The explicit expressions for E^tt and E^tt are presented in (<ref>) and (<ref>).






§ SOLUTIONS IN STAROBINSKY'S MODEL




Even the reduced system of equations (<ref>), (<ref>) is too complicated to be solved analytically.
An obvious strategy is to make additional assumptions, such as e.g. a constant scalar curvature R=R_0, a particular relation between B and A, or a specific ansatz for the functions B and A. 
While such assumptions lead to strong simplifications and, in some cases, admit exact analytical solutions, they do not provide a general strategy to systematically explore the space of solutions.
In this letter, we do not make any additional simplifying assumptions, but instead classify different solutions according to their asymptotic behaviour at small and large radii and then connect these two regimes by solving the system (<ref>), (<ref>) numerically. 
Various special radii r, corresponding to relevant regimes in the problem, are listed in Table <ref>.







 §.§ Asymptotics at the origin




Following the strategy of <cit.>, applied in the context of QG in <cit.>, the solutions can be systematically classified by expanding B and A in a Frobenius series around r=0,

    B(r)=   b_t(r^t+b_t+1r^t+1+…),
    
    A(r)=   a_sr^s+a_s+1r^s+1+… .
 
In (<ref>), b_t has been factored out, as it can be absorbed by rescaling the time coordinate t→t̃=√(b_t)t.
Consistent combinations of (s,t) follow from the indicial equations, which are derived by inserting (<ref>), (<ref>) into (<ref>), (<ref>).
We find three subcases:

    s<0:     no solution,
    
    s=0:     t arbitrary,
    
    s>0:     s=(t^2+2t+4)(t+4)^-1.
 
For s<0, no solution exists.
For s=0, the value of t is unconstrained and there are infinitely many solutions with fixed a_0=(1+t/2+t^2/4).
For s>0 and s,t∈ℝ, there are infinitely many solutions of (<ref>).
For s>0 and s,t∈ℤ, there are only six different integer families (s,t) which satisfy (<ref>).[In QG, only the (0,0), (-1,1) and (2,2) families were found <cit.>.]
The numerical analysis shows that, in contrast to the s>0 families, the s=0 families cannot be connected to asymptotically flat solutions of the full non-linear equations, except for t=0.
Therefore, there are a total of seven asymptotically flat, integer Frobenius families (s,t) summarized in Table <ref>.
In the limit r→0 for s≠0, the leading r-dependence of the Kretschmann scalar K:=R_μνρσR^μνρσ∝ r^-2(2+s) implies that theses solutions have a singular geometry.[The (0,0) family is special.
Extracting the Frobenius coefficients in (<ref>), (<ref>) up to 𝒪(r^8), leads to two solutions r→0limK= (2M_S^4/3) [4+ã_2(1+ã_2)∓2(4+ã_2(2+ã_2))^1/2] with ã_2:=6a_2/M_S^2 a free parameter. 
For ã_2=0, the minus sign (corresponding to Minkowski space) implies r→0limK=0. 
]




The vacuum solutions of GR are also solutions of (<ref>).
In particular, Minkowski space (MS) is included in the (0,0) family, while the Schwarzschild solution (SS) is included in the (1,-1) family.





 §.§ Linearised solutions and asymptotic flatness




The linearised EOMs in the weak field regime r≥ r_L are derived by inserting the decomposition

    B(r)=1+V(r),     A(r)=1+W(r),

into (<ref>),(<ref>) and by keeping only terms linear in V or W.
The linearised EOMs have the analytic solution <cit.>,

    V(r)=   C+C_2/r+∑_σ=-,+C_0^σe^σ M_Sr/r, 
    
    W(r)=   -C_2/r+∑_σ=-,+C_0^σe^σ M_Sr/r(1-σ M_Sr).

The four integration constants C, C_2, and C_0^± agree with the result D=4 found 
by the DAE analysis (<ref>).[The same DAE analyses in QG and EWG lead to D=6 and D=4, matching the number of integration constants of the linearised equations, i.e. twice the number of propagating particles.]
Asymptotic flatness of the geometry requires

    lim_r→∞B(r)=1,    lim_r→∞A(r)=1.

The limits (<ref>), in turn, require the constant C and the constant C_0^+ of the rising scalaron Yukawa potential to vanish, reducing the number of free parameters to two, 

    V^∞(r):=   C_2/r+C_0^-e^-M_Sr/r,
    
    W^∞(r):=   -C_2/r+C_0^-e^-M_Sr/r(1+M_Sr).

Matching (<ref>) with the Newtonian potential at r→∞, 

    lim_r→∞V^∞(r)=-2G_NM_T^∞/r,

relates C_2 with the total mass M_T^∞ at infinity.
In terms of Planck units, this implies C_2=-16π M_T^∞.
In contrast, C_0^- cannot be directly related to an observable at infinity and a priori remains an undetermined parameter.





 §.§ Numerical algorithm




We connect the asymptotics of the solutions A and B in the Frobenius regime with those in the linearised regime, by numerically solving the full non-linear equations (<ref>).
This is a boundary problem which is difficult to solve – even numerically.
Therefore, we use a double-shooting algorithm.
We first formulate the problem as initial value problem (IVP) with “asymptotic initial conditions” generated by matching the full numerical solution to the linearised solution at large radii r_I≥ r_L, 

    B(r_I)=   1+V^∞(r_I),     B_1(r)|_r_I=V^∞_1(r)|_r_I,
      A(r_I)=   1+W^∞(r_I),     A_1(r)|_r_I=W^∞_1(r)|_r_I.

For fixed M_S, the IVP only depends on the two parameters C_0^- and C_2.
Fixing C_0^- and C_2, we integrate the system inwards from the initial radius r_I to the final radius r_E at which the solutions A and B are dominated by the leading term of the Frobenius series.
Matching the numerical solutions A and B at 0<r_E≤ r_F to (<ref>) and (<ref>), the numerical values (s_N,t_N) are extracted from the slope of a linear log-log fit with constants c_B and c_A,

    ln B(r_E) =   c_B+ t_Nln r_E+…,
    ln A(r_E)  =   c_A+ s_Nln r_E +…,

In this way, for each pair of initial values (C_2,C_0^-), a pair of numerical values (s_N,t_N) is found.
Repeating this procedure by scanning over the (C_2,C_0^-) parameter space, we only store those values of (C_2,C_0^-) for which the numerical values (s_N,t_N) coincide with the (s,t) values of one of the integer Frobenius families up to a tolerance max{Δ s,Δ t}< δ_S=10^-3, with the fractional differences  Δ s:=(s-s_N)/s and Δ t:=(t-t_N)/t, between the corresponding numerical values (s_N,t_N) and the integer Frobenius values (s,t).

By construction of the shooting algorithm, the numerical solutions A and B are only strictly valid up to r_I.
Despite the fact that M_S parametrically enters the EOMs (<ref>),(<ref>), from the point of view of the initial conditions (<ref>),(<ref>), the value of r_I is tied to the value of M_S, which we choose according to the CMB constraints (<ref>), i.e. in Planck units M_S=10^-5.
For a given r_I, different “confidence regions” in (C_2,C_0^-) space correspond to different tolerance values max{|V^∞(r_I)|,|W^∞(r_I)|}<δ_L, which quantify how well the linearised approximation is satisfied at r_I.
While a large r_I would be desirable as it allows for larger values of (C_2,C_0^-) within the same confidence interval, it is computationally expensive.
For the numerical analysis, we choose r_I=10 (in Planck units).
Likewise, choosing r_E too small, the singularity of the EOMs at r=0 might cause the solution to diverge already for r>r_E, while choosing r_E too large might result in a poor fit (<ref>), (<ref>).
We choose the default value r_E=10^-4, but if the quality of the fit falls below a tolerance δ_F, the code dynamically adjusts by automatically lowering the value of r_E. 

We performed various consistency checks.
First, we repeated the Frobenius analysis also in Kundt coordinates (as done in <cit.> for QG) and found the same integer Frobenius families as in SS coordinates.
Second, as a check of our numerical implementation, we applied our algorithm to the model of QG and found good agreement with the numerical results obtained in <cit.>.
Third, we checked that the numerical solutions of the reduced system (<ref>),(<ref>) also satisfy the original system (<ref>), (<ref>), by re-substituting them into the original system.
Fourth, we checked that for all (including non-integer) numerically found values (s_N,t_N), the Frobenius constraint (<ref>) is satisfied with a tolerance not exceeding |s-(t^2+2t+4)(t+4)^-1|<δ_C=10^-3.[The numerically extracted (s,t) values, all satisfying (<ref>), vary smoothly under variation of (C_2,C_0^-), implying that non-integer families are qualitatively not different from integer ones.]





§ RESULTS




We find that all integer Frobenius families can be connected to asymptotically flat solutions and lead to the “phase-diagram” Fig. <ref>.


Except for the (0,0) family, which appears as a single point at the origin in Fig. <ref>, all integer Frobenius families trace out lines in parameters space, which emanate in the vicinity of the origin and expand outwards to larger (C_2,C_0^-) values without intersecting. While Minkowski space is contained in the (0,0) family, the Schwarzschild solution, which is formally part of the (1,-1) family, is not contained in the numerically found solutions because it has a horizon.[Numerically, we could neither find solution with s<0, confirming the expected result from (<ref>), nor could we find any of the (0,t) except for the (0,0) family, suggesting that the (0,t), t≠0 families cannot be connected to asymptotically flat solutions.
The (0,0) family appears as a single point in Fig. <ref>, as to all orders (<ref>), (<ref>) depends on a single Frobenius parameter.]

According to Table <ref>, all integer Frobenius families (except the (0,0) family) have a singular geometry.
By construction, none of the found solution has a horizon, such that the corresponding solutions all feature a naked singularity at the origin.[In the shooting algorithm solutions with horizons are discarded as they would diverge at some r_H>r_E such that the numerical integration would fail before a successful fit of (s_N,t_N) could be obtained.
In any case, the uniqueness theorem for Starobinsky's model guarantees that the SS is the only spherically symmetric, asymptotically flat vacuum solution with a horizon <cit.>.
Nevertheless, by adapting the shooting algorithm, making it capable of detecting a horizon, we numerically “confirm” this, as, similar to the findings in <cit.>, the only positive mass solution with horizon we found is the SS solution lying on the negative C_2 axis in Fig. <ref>.
]

Since, according to (<ref>), C_2 is connected to the total mass at infinity, positivity of M_T^∞ requires C_2 to be negative.
The physical significance of the negative mass families (1,-1), (1,0) and (2,-2) remains unclear, but we focus on the positive mass families (2,2), (7,-3) and (7,8) which lie in the negative C_2 half-plane of Fig. <ref>. 


In contrast to C_2, the parameter C_0^- cannot directly be related to an observable at infinity.
However, the phase diagram Fig. <ref> shows that for each integer Frobenius family, C_0^- can uniquely be expressed as a function of C_2 and therefore as a function of M_T^∞. 
Within the reliable C_2 intervals of each family, the function C_0^-(C_2) is well described by a quadratic fit shown in Fig. <ref> for the positive mass families.
This implies that all solutions in the integer Frobenius families are completely characterized by two (s,t) values and by one single physical parameter M_T^∞.

In the numerical analysis we assumed M_S=10^-5 and r_I=10, such that, even for radii r≥ r_L well within the linearised regime, we can expand the exponential in (<ref>),

    V^∞(r)=C_2+C_0^-/r+O(M_Sr),     r≤ r_S,
 
which implies the constraint (C_2+C_0^-)/r_I<δ_L.
Using the quadratic fits obtained in Fig. <ref> with the identification C_2=-16 M_T^∞, we encounter the interesting situation that even in the linearised regime, the “effective quasi-local mass” C_2+C_0^- might vanish or become negative, rendering the weak field potential V^∞(r) repulsive.
For the (2,2) family, we obtain in this way

    V^∞(r)≈-16π M_T^∞/r(1.76-0.56 M_T^∞),
 
which has a zero at M_T^∞≈π and even becomes negative for larger values of M_T^∞.
This scalaron-induced “mass screening mechanism” for r<r_S might have interesting implications for the  solutions, as e.g. avoiding a repulsive potential would imply a constraint on the total mass at infinity M_T^∞<π, which, however, in the present case for M_T^∞=π lies outside the region of validity of the quadratic fit of the (2,2) family as -16π^2=C_2≪-1.2, cf. Fig. <ref>.
Similar considerations for the (7,-3) and (7,8) families would imply an anti-screening. 

Summarizing, we found new static, spherically symmetric and asymptotically flat vacuum solutions of Starobinsky's quadratic f(R) model which feature naked singularities.
The phase diagram Fig. <ref> might be interpreted as  “Frobenius no-hair theorem” as all solutions in the seven integer Frobenius families labelled by the two Frobenius indices (s,t) are characterized by one single physical parameter – the total mass at infinity M_T^∞.
Besides the condition of positive total mass at infinity, we found a mass (anti)screening mechanism, which depends on the Frobenius family and the scalaron mass M_S, and might lead to additional constraints on the Frobenius solutions.  
The implications of the naked singularity solutions in the context of the black hole information paradox requires a more detailed study involving the analysis of the shape of different solutions within a given Frobenius family. We hope to address this question in a forthcoming work.
Finally, it would be interesting to apply the presented methods to other relevant f(R) models and modifications of GR.





§ EXPLICIT FORM OF EOMS



The explicit expressions for the coefficients ℰ^tt and ℰ^rr of (<ref>), which enter the EOMs (<ref>) and (<ref>), read

    ℰ^tt=   
    M_P^2/96 M_S^2 r^4
    	A^5 B^3{A^3 [16 B^4 (3 M_S^2 r^3 A_1+5)-4 r^2 B^2 (4
    	r B_1 {11 B_2+3 r B_3}+4 B_1^2+9 r^2 B_2^2)..
    
    	   .+4 r^3 B
    	B_1^2 (26 B_1+29 r B_2)-49 r^4 B_1^4+16 r^3 (4 B_3+r
    	B_4) B^3]-2 r A^2 B [4 r B^2 (r {[6 A_2+r
    	A_3] B_1...
    
    	   ..+4 r A_2 B_2}-2 A_1 {B_1-r (8 B_2+3 r
    	B_3)})-6 r^2 B B_1 (2 B_1 {5 A_1+r A_2}+9 r
    	A_1 B_2)+29 r^3 A_1 B_1^3
    
    	   .+16 B^3 (r {A_2+r A_3}-2
    	A_1)]+r^2 A B^2 A_1 [4 r B (B_1 {28 A_1+13 r
    	A_2}+19 r A_1 B_2).
    
    	   .-57 r^2 A_1 B_1^2+16 B^2 (5 A_1+13 r
    	A_2)]-56 r^3 B^3 A_1^3 [r B_1+4 B]+16 A^5 B^4 [3
    	M_S^2 r^2+1]
    
    	   .-48 A^4 B^4 [M_S^2 r^2+2]},
    
    
    ℰ^rr=   
    M_P^2/96 M_S^2 r^4 A^3 B^4{2 r^2 A B [r B_1+4 B] [2 B (B_1 (4
    	A_1+r A_2)+2 r A_1 B_2)-3 r A_1 B_1^2+8 B^2 A_2].
    
    	   -7 r^2
    	B^2 A_1^2 [r B_1+4 B]^2+48 A^3 B^3 [M_S^2 r^3 B_1+B
    	(M_S^2 r^2-2)]+A^2 [-32 r B^3 (r {2 B_2+r
    	B_3}-2 B_1).
    
    	   .+4 r^2 B^2 (-2 r B_1 {r B_3-6
    	B_2}+16 B_1^2+r^2 B_2^2)+4 r^3 B B_1^2 (3 r B_2-4
    	B_1)-7 r^4 B_1^4+112 B^4]
    
    	   .-16 A^4 B^4 [3 M_S^2
    	r^2+1]}.

The explicit expressions for the coefficients E^tt and E^rr, which enter the EOMs (<ref>) and (<ref>), read

    E^tt=   
    M_P^2/8 M_S^2 r^4 A^4 B^2 {r B_1+4
    	B}^2
    
    	   ×{A^2 [-4 r^2 B^3 (B_1^2 {M_S^2 r^3 A_1+4}+8 r
    	B_1 B_2+2 r^2 B_2^2)+8 r B^4 B_1 (r {2-M_S^2 r^2}
    	A_1+2)+32 B^5 (r {M_S^2 r^2....
    
    	   ...+2} A_1+1)-4 B^2 (r^5
    	B_1 B_2^2-4 r^3 B_1^3)+2 r^4 B B_1^3 (B_1+2 r B_2)-r^5
    	B_1^5]-2 r A B A_1 [r B_1+4 B] [-2 r^2 B B_1
    	(B_1..
    
    	   ..+r B_2)+r^3 B_1^3+8 B^3]+r^2 B^2 A_1^2 [2 B-r
    	B_1] [r B_1+4 B]^2-4 A^3 B^2 [2 r^2 B B_1 ({3
    	M_S^2 r^2+2} B_1-M_S^2 r^3 B_2).
    
    	   .-8 B^2 (r {M_S^2
    	r^2-1} B_1+M_S^2 r^4 B_2)+M_S^2 r^5 B_1^3+8 B^3 (M_S^2
    	r^2+2)]+16 A^4 B^4 [r (2 M_S^2 r^2+1) B_1.
    
    	   ..+2 B
    	(M_S^2 r^2+1)]},
    
    
    E^rr=   
    -M_P^2/96
    	M_S^2 A^3 B^4 r^5 {2 M_S^2 A^2 r [4 B+r B_1] B^2+r
    	A_1 B_1 [4 B+r B_1] B+A [r^2 B_1^3-8 B^2 B_1-2 B
    	r (2 B+r B_1) B_2]}
    
    	   ×{32 M_S^2 A^6 r^2 [3 M_S^2 r^2+1] [4 B+r
    	B_1] B^6+16 A^5 [8 (-3 M_S^4 r^4+10 M_S^2 r^2+8) B^3+2 r
    	({-15 M_S^4 r^4+2 M_S^2 r^2....
    
    	   ....+12} B_1-2 r {5 M_S^2 r^2+1}
    	B_2) B^2-2 r^2 B_1 (M_S^2 B_2 r^3+{3 M_S^4 r^4+1}
    	B_1) B-r^3 (5 M_S^2 r^2+3) B_1^3] B^4.
    
    	   .+r^3 A_1^3 [16
    	B-r B_1] [4 B+r B_1]^3 B^3-A r^2 A_1 [4 B+r
    	B_1]^2 [32 (2 A_1+r A_2) B^3+4 r (2 B_1 {9 A_1+r
    	A_2}...
    
    	   ..+5 r A_1 B_2) B^2-2 r^2 A_1 B_1 (B_1+2 r B_2)
    	B-r^3 A_1 B_1^3] B^2-2 A^4 [-7 M_S^2 B_1^5 r^7+8 B B_1^3
    	(2 M_S^2 r^3 B_2..
    
    	   .-{M_S^2 r^2+2} B_1) r^4+4 B^2 B_1
    	(-M_S^2 B_2^2 r^4+6 M_S^2 B_1 B_2 r^3+{54 M_S^2 r^2+4}
    	B_1^2) r^3+8 B^3 (-2 M_S^2 B_2^2 r^4.
    
    	   .-34 M_S^2 B_1 B_2
    	r^3-{8 M_S^2 r^2+[3 M_S^2 r^3+r] A_1-12} B_1^2) r^2+16
    	B^4 ({3 M_S^2 r^2-2 (3 M_S^2 r^3+r) A_1+56} B_1.
    
    	   ..+2 r
    	{6-5 M_S^2 r^2} B_2) r+64 B^5 (11 M_S^2 r^2+16)]
    	B^2+A^2 r [4 B+r B_1] [32 (A_1 {r (11 M_S^2
    	r^2+8) A_1-8}..
    
    	   .+4 r A_2) B^5+16 r (B_1 {A_1 [r (5
    	M_S^2 r^2+4) A_1-3]+10 r A_2}+2 r {2 r A_2-A_1}
    	B_2) B^4-2 (16 A_1 B_2^2 r^4.
    
    	   .-8 B_1 A_2 B_2 r^4+A_1
    	.{M_S^2 r^3 A_1-48} B_1^2 r^2) B^3-4 r^4 B_1 (2 A_2
    	B_1^2-10 A_1 B_2 B_1+r A_1 B_2^2) B^2-4 r^4 A_1 B_1^3 (7
    	B_1.
    
    	   ..+2 r B_2) B+5 r^5 A_1 B_1^5] B+A^3 [-256 ({3
    	M_S^2 r^2+2} A_2 r^2-4 A_1 r-4) B^7+64 r (B_1 {r
    	[(6-9 M_S^2 r^2) A_1....
    
    	   ...-2 r (3 M_S^2 r^2+2)
    	A_2]+22}+r {r (3 M_S^2 r^2+2) A_1+7} B_2)
    	B^6+16 r^2 B_1 (r^2 {3 M_S^2 r^2+2} A_1 B_2.
    
    	   .-B_1 {r
    	[(39 M_S^2 r^2+22) A_1+r (3 M_S^2 r^2+2)
    	A_2]-14}) B^5-8 r^3 ({3 r [3 M_S^2 r^2+4]
    	A_1-26} B_1^3+8 r B_2 B_1^2.
    
    	   .+4 r^2 B_2^2 B_1-2 r^3 B_2^3) B^4+4
    	r^4 B_1^2 (3 {M_S^2 A_1 r^3+16} B_1^2+40 r B_2 B_1-2 r^2
    	B_2^2) B^3+4 r^5 B_1^3 (-14 B_1^2+7 r B_2 B_1.
    
    	   ...+3 r^2 B_2^2)
    	B^2-2 r^6 B_1^5 (7 B_1+6 r B_2) B+3 r^7 B_1^7]}.






















