# 2nd order perturbations applied to f(\chi)=\chi^{3/2}

Copyright (C) 2013 T. Bernal <tbernal@astro.unam.mx>, S. Mendoza <sergio@mendozza.org>
and L.A. Torres <ltorres@astro.unam.mx>

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA

For information on GNU (or the Free Software Foundation -FSF) visit
http://www.gnu.org.  For information on the GNU General Public Licence,
visit: http://www.gnu.org/licences/licences.html#GPL

We first reset and kill(all) variables.  This is to start a new session:

reset();
kill(all);

First load the ctensor package with the csetup() interface. That asks us the way
we want to input the metric we would like to build. 
In order of appearance, the answers are:
- 4 dimensions
- y in order to change the coordinate names (see following line)
- r,theta,phi,t coordinates (Note: time coordinate t is the 4th coordinate)
- y to introduce a new metric
- y for a diagonal metric
- g_11 is introduced
- g_22 is introduced
- g_33 is introduced
- g_44 is introduced (Note that this is our g_00)
- We assign dependencies of the metric coefficients as functions of the coordinates.
- y to display the metric g_{\mu \nu} in matrix form
Note that in the following epsilon = 1 / c, where c is the speed of light (we are
going to make perturbations and this definition is coherent).

load(ctensor);
csetup();
4;
y;
[r,theta,phi,t];
1;
1;
- 1 + epsilon^2 * g112;
- r^2;
- r^2 * sin(theta)^2;
1 + epsilon^2 * g002 + epsilon^4 * g004;
depends([g112,g002,g004],r);
y;
print(Continue);

We now define by hand the g^{\mu\nu}'s.

gup[4,4] : 1 - epsilon^2 * g002 - epsilon^4 * g004;
gup[1,1] : -1 - epsilon^2 * g112;
gup[2,2] : ug[2,2];
gup[3,3] : ug[3,3];

Show Christoffel symbols (note that first two indices are the symmetric ones):

christof(mcs);

We calculate the Christoffel symbols at zeroth order O(0) by taking the limit when the
velocity of light c tends to infinity, which is the same as taking the limit when
epsilon tends to zero. This is the value epsilon one gets in spherical coordinates in
flatland.

a1 : limit(mcs[1,1,1],epsilon,0);
a2 : limit(mcs[1,2,2],epsilon,0);
a3 : limit(mcs[1,3,3],epsilon,0);
a4 : limit(mcs[1,4,4],epsilon,0);
a5 : limit(mcs[2,2,1],epsilon,0);
a6 : limit(mcs[2,3,3],epsilon,0);
a7 : limit(mcs[3,3,1],epsilon,0);
a8 : limit(mcs[3,3,2],epsilon,0);
a9 : limit(mcs[4,4,1],epsilon,0);

The 2nd order O(2) of the Christoffel symbols is now calculated by using the results
of the Christoffel symbols at O(0). This is done in two steps. The first step
calculates just the 2nd order (or equivalently, the 2nd order minus the zeroth order):

b1 : limit(epsilon^(-2) * (mcs[1,1,1] - a1), epsilon,0);
b2 : limit(epsilon^(-2) * (mcs[1,2,2] - a2), epsilon,0);
b3 : limit(epsilon^(-2) * (mcs[1,3,3] - a3), epsilon,0);
b4 : limit(epsilon^(-2) * (mcs[1,4,4] - a4), epsilon,0);
b5 : limit(epsilon^(-2) * (mcs[2,2,1] - a5), epsilon,0);
b6 : limit(epsilon^(-2) * (mcs[2,3,3] - a6), epsilon,0);
b7 : limit(epsilon^(-2) * (mcs[3,3,1] - a7), epsilon,0);
b8 : limit(epsilon^(-2) * (mcs[3,3,2] - a8), epsilon,0);
b9 : limit(epsilon^(-2) * (mcs[4,4,1] - a9), epsilon,0);

The 4th order O(4) of the Christoffel symbols is now calculated by using the results
of the Christoffel symbols at O(0) and O(2). This step calculates just the 4th order
(or equivalently, the 4th order minus the 2nd and the zeroth order):

c1 : limit(epsilon^(-4) * fullratsimp(mcs[1,1,1] - epsilon^2*b1 - a1), epsilon,0);
c2 : limit(epsilon^(-4) * fullratsimp(mcs[1,2,2] - epsilon^2*b2 - a2), epsilon,0);
c3 : limit(epsilon^(-4) * fullratsimp(mcs[1,3,3] - epsilon^2*b3 - a3), epsilon,0);
c4 : limit(epsilon^(-4) * fullratsimp(mcs[1,4,4] - epsilon^2*b4 - a4), epsilon,0);
c5 : limit(epsilon^(-4) * fullratsimp(mcs[2,2,1] - epsilon^2*b5 - a5), epsilon,0);
c6 : limit(epsilon^(-4) * fullratsimp(mcs[2,3,3] - epsilon^2*b6 - a6), epsilon,0);
c7 : limit(epsilon^(-4) * fullratsimp(mcs[3,3,1] - epsilon^2*b7 - a7), epsilon,0);
c8 : limit(epsilon^(-4) * fullratsimp(mcs[3,3,2] - epsilon^2*b8 - a8), epsilon,0);
c9 : limit(epsilon^(-4) * fullratsimp(mcs[4,4,1] - epsilon^2*b9 - a9), epsilon,0);

Now we write the Christoffel symbols up to the fourth order:

mcs111 : a1 + epsilon^2 * b1 + epsilon^4 * c1;
mcs122 : a2 + epsilon^2 * b2 + epsilon^4 * c2;
mcs133 : a3 + epsilon^2 * b3 + epsilon^4 * c3;
mcs144 : a4 + epsilon^2 * b4 + epsilon^4 * c4;
mcs221 : a5 + epsilon^2 * b5 + epsilon^4 * c5;
mcs233 : a6 + epsilon^2 * b6 + epsilon^4 * c6;
mcs331 : a7 + epsilon^2 * b7 + epsilon^4 * c7;
mcs332 : a8 + epsilon^2 * b8 + epsilon^4 * c8;
mcs441 : a9 + epsilon^2 * b9 + epsilon^4 * c9;

The determinant of the metric tensor g is given by:

g:gdet;

We calculate Ln(sqrt(-|g|)) up to the 4th order:

difflogg : diff(log(sqrt(-g)),r);
dlogg0 : limit(difflogg,epsilon,0);
dlogg2 : limit(epsilon^(-2) * (difflogg - dlogg0), epsilon, 0);
dlogg4 : limit(epsilon^(-4) * fullratsimp(difflogg - epsilon^2*dlogg2 - dlogg0), epsilon, 0);
dlogg : dlogg0 + epsilon^2 * dlogg2 + epsilon^4 * dlogg4;

We define the dependence of the Ricci scalar up to the 4th order as function of the
radius, and the dependence of the Ricci tensor components R00 and R22 up to the 4th
order as functions of the radius:

depends([R2,R4,R002,R222,R004,R224],r);

We define the Ricci scalar up to the 4th order, the function f(R) we want to use and
its derivatives:

R : epsilon^2 * R2 + epsilon^4 * R4 ;
f   : epsilon^(2*b) * R2^b;
fp  : b * epsilon^(2*(b-1)) * R2^(b-1) ;
f2p : b * (b-1) * ( epsilon^(2*(b-2))*R2^(b-2) + (b-2)*epsilon^(2*(b-3))*R2^(b-3)*(epsilon^4*R4) ) ;
f3p : b * (b-1) * (b-2) * ( epsilon^(2*(b-3))*R2^(b-3) + (b-3)*epsilon^(2*(b-4))*R2^(b-4)*(epsilon^4*R4) );

We define the following factors, A and B, to calculate the trace equation:

Adef : (diff(gup[1,1],r) + gup[1,1] * dlogg) * diff(R,r) + gup[1,1] * diff(R,r,2) ;
Bdef : expand(gup[1,1]*(diff(R,r))^2) ;
A2 : limit(epsilon^(-2) * Adef, epsilon, 0) ;
A4 : limit(epsilon^(-4) * fullratsimp(Adef - epsilon^2*A2), epsilon, 0) ;
A6 : limit(epsilon^(-6) * fullratsimp(Adef - epsilon^2*A2 - epsilon^4*A4), epsilon, 0) ;
B4 : limit(epsilon^(-4) * Bdef, epsilon, 0) ;
B6 : limit(epsilon^(-6) * fullratsimp(Bdef - epsilon^4*B4), epsilon, 0) ;
B8 : limit(epsilon^(-8) * fullratsimp(Bdef - epsilon^4*B4 - epsilon^6*B6), epsilon, 0) ;
A : epsilon^2*A2 + epsilon^4*A4 + epsilon^6*A6;
B : epsilon^4*B4 + epsilon^6*B6 + epsilon^8*A8;

We define the factor H to calculate the trace equation:

H : factor(fullratsimp( 3 * f2p * A + 3 * f3p * B )) ;

Trace equation (in vacuum):

G : fp * R - 2 * f + H ;

Trace equation at the order (2b-2):

G1 : factor(limit(epsilon^(-(2*b-2))*G, epsilon, 0)) ;
g1 : expand(G1 / (-3*b*(b-1)*R2^(b-3))) ;

Trace equation at the order (2b):

g3 : limit(epsilon^(-2*b) * fullratsimp(G - epsilon^(2*b-2)*G1), epsilon, 0);

Now we make some algebra to reduce the expression:

l3 : expand(fullratsimp(factor(g3))) ;
m3 : expand(l3 / R2^(b-4));
n3 : expand(b3 + 3*b*(b-1) * ((b-2)*R4*g1 + g112*R2*g1) );

We define the following factors to obtain the 00 component of the field equations:

Cdef : diff(R,t,2) - mcs441 * diff(R,r) - lg[4,4] * ( (diff(gup[1,1],r)
    + gup[1,1] * dlogg) * diff(R,r) + gup[1,1] * diff(R,r,2) ) ;
Ddef : (diff(R,t))^2 - lg[4,4] * gup[1,1] * (diff(R,r))^2 ;
C2 : limit(epsilon^(-2) * fullratsimp(Cdef), epsilon, 0) ;
C4 : limit(epsilon^(-4) * fullratsimp(Cdef - epsilon^2*C2), epsilon, 0) ;
C6 : limit(epsilon^(-6) * fullratsimp(Cdef - epsilon^2*C2 - epsilon^4*C4), epsilon, 0) ;
D4 : limit(epsilon^(-4) * Ddef, epsilon, 0) ;
D6 : limit(epsilon^(-6) * fullratsimp(Ddef - epsilon^4*D4), epsilon, 0) ;
D8 : limit(epsilon^(-8) * fullratsimp(Ddef - epsilon^4*D4 - epsilon^6*D6), epsilon, 0) ;
C : epsilon^2*C2 + epsilon^4*C4 + epsilon^6*C6 ;
D : epsilon^4*D4 + epsilon^6*D6 + epsilon^8*D8 ;

We define the factor H00 to obtain the field equation G00:

H00 : factor(fullratsimp(- f2p * C - f3p * D)) ;

Field equation 00 (in vacuum):

G00 : fp * (epsilon^2*R002) - f * lg[4,4] / 2 + H00 ;

Field equation 00 at the order (2b-2):

G001 : factor(limit(epsilon^(-(2*b-2))*G00, epsilon, 0));

Field equation 00 at the order (2b):

G003 : limit(epsilon^(-2*b) * fullratsimp(G00 - epsilon^(2*b-2)*G001), epsilon, 0);
j: expand(G003 * R2^(-b+4));
k: expand(j + b*(b-1)*R2*(g112+g002)*g1 + b*(b-1)*(b-2)*R4*g1) ;
l : expand( (k - n3/3) / R2^2 );

We define the following factors to obtain the theta-theta component of the field
equations:

Edef : diff(R,theta,2) - mcs221 * diff(R,r) - lg[2,2] * ( (diff(gup[1,1],r)
    + gup[1,1] * dlogg) * diff(R,r) + gup[1,1] * diff(R,r,2) ) ;
Fdef : (diff(R,theta))^2 - lg[2,2] * gup[1,1] * (diff(R,r))^2 ;
E2 : limit(epsilon^(-2) * fullratsimp(Edef), epsilon, 0) ;
E4 : limit(epsilon^(-4) * fullratsimp(Edef - epsilon^2*E2), epsilon, 0) ;
E6 : limit(epsilon^(-6) * fullratsimp(Edef - epsilon^2*E2 - epsilon^4*E4), epsilon, 0) ;
F4 : limit(epsilon^(-4) * Fdef, epsilon, 0) ;
F6 : limit(epsilon^(-6) * fullratsimp(Fdef - epsilon^4*F4), epsilon, 0) ;
F8 : limit(epsilon^(-8) * fullratsimp(Fdef - epsilon^4*F4 - epsilon^6*F6), epsilon, 0) ;
E : epsilon^2*E2 + epsilon^4*E4 + epsilon^6*E6 ;
F : epsilon^4*F4 + epsilon^6*F6 + epsilon^8*F8 ;

We define the factor H22 to obtain the field equation G22:

H22 : factor(fullratsimp(- f2p * E - f3p * F)) ;

Field equation 22 (in vacuum):

G22 : fp * (epsilon^2*R222) - f * lg[2,2] / 2 + H22 ;

Field equation 22 at the order (2b-2):

G221 : factor(limit(epsilon^(-(2*b-2))*G22, epsilon, 0));

Field equation 22 at the order (2b):

G223 : limit(epsilon^(-2*b) * fullratsimp(G22 - epsilon^(2*b-2)*G221), epsilon, 0);
n: expand( G223 *3 / (R2^(b-4)*r^2) );
o: expand( n/3 + n3/3 - g112*b*(b-1)*R2*g1 - b*(b-1)*(b-2)*R4*g1 );

The equation for R2 can be written as follows:

'diff(R2,r) * R2^(b-2) = Ct/r^2 ;
RR2 : ode2(%,R2,r);
nonzero;
print(Continue);

And the solution for the Ricci scalar at the second order is given by:

solve(RR2,R2);
yes;
print(Continue);

The solution for the Ricci scalar at the second order for b=3/2 and very far from
the source mass (%c*Ct=0) is R2=Rt/r^2, for Rt=const.