> restart:

> grtw():

> read `/usr/local/MapleVR5/Grtii(5)/Lib/pertutils.mpl`;

> qload(qschw):

As can be seen from the form of the metric, the lower case functions are the first order perturbations, and the upper case functions are the second order perturbations. To begin our analysis we again calculate the exact contravariant metric tensor, which is then quadratically perturbed using `` quadpert() ''.

> grcalc(g(up,up));

> gralter(g(up,up),quadpert,simplify,factor);

Component simplification of a GRTensorII object:

Applying routine quadpert to object g(up,up)

Applying routine simplify to object g(up,up)

Applying routine factor to object g(up,up)

> grdisplay(g(up,up));

With these quadratically perturbed contra/co-variant components of the metric we can now calculate the perturbed Ricci tensor to second order.

> grcalcalter(R(dn,dn),13);

Simplification will be applied during calculation.

Applying routine `Apply constraints repeatedly` to object g(dn,dn,pdn)

Applying routine `Apply constraints repeatedly` to object Chr(dn,dn,dn)

Applying routine `Apply constraints repeatedly` to object Chr(dn,dn,up)

Applying routine `Apply constraints repeatedly` to object R(dn,dn)

> gralter(R(dn,dn),expand);

Component simplification of a GRTensorII object:

Applying routine expand to object R(dn,dn)

> grdisplay(R(dn,dn));

For completeness, the definitions of the Legendre function, P_2(cos(theta)) , that we will need to calculate the projections listed in Table 1, and the standard Zerilli substitutions, that we will need to eliminate H_1 , and, Kdot , appear below.

> p2 := 5/2*(3*cos(theta)^2-1)/2;

> newH(r,t) := (2*r^2-6*m*r-3*m^2)/(r-2*m)/(2*r+3*m)*chi(r,t)+r^2/(r-2*m)*eta(r,t);

> Kdot(r,t) := 6*(r^2+m*r+m^2)/r^2/(2*r+3*m)*chi(r,t)+eta(r,t);

Following the prescription set out in Table 1, we can now extract the linear contributions to the seven non-trivial Ricci components.

> a1:=coeff(collect(grcomponent(R(dn,dn),[theta,theta]),epsilon),epsilon,1):

> g1:=coeff(collect(grcomponent(R(dn,dn),[phi,phi]),epsilon),epsilon,1):

> ag1:=int(expand(a1-g1/sin(theta)^2)*(2*cot(theta)*diff(p2,theta)+2*(2+1)*p2)*sin(theta),
theta=0..Pi);

> leqn[1]:= {h[2](r,t) = solve(ag1,h[2](r,t))};

> a2:=subs(leqn[1],a1):

> g2:=subs(leqn[1],g1):

> grmap(R(dn,dn),subs,leqn[1],'x');

Applying routine subs to R(dn,dn)

> gralter(R(dn,dn),simplify,expand);

Component simplification of a GRTensorII object:

Applying routine simplify to object R(dn,dn)

Applying routine expand to object R(dn,dn)

> b1:=coeff(collect(grcomponent(R(dn,dn),[theta,r]),epsilon),epsilon,1):

> leqn[2]:=termsimp(collect(1/3*int(b1*diff(p2,theta)*sin(theta),theta=0..Pi),lperts));

> c1:=coeff(collect(grcomponent(R(dn,dn),[theta,t]),epsilon),epsilon,1):

> leqn[3]:=termsimp(collect(1/3*int(c1*diff(p2,theta)*sin(theta),theta=0..Pi),lperts));

> d1:=coeff(collect(grcomponent(R(dn,dn),[t,r]),epsilon),epsilon,1):

> leqn[4]:=termsimp(collect(int(d1*p2*sin(theta),theta=0..Pi),lperts));

> e1:=coeff(collect(grcomponent(R(dn,dn),[r,r]),epsilon),epsilon,1):

> leqn[5]:=termsimp(collect(int(e1*p2*sin(theta),theta=0..Pi),lperts));

> f1:=coeff(collect(grcomponent(R(dn,dn),[t,t]),epsilon),epsilon,1):

> leqn[6]:=termsimp(collect(int(f1*p2*sin(theta),theta=0..Pi),lperts));

> leqn[7]:=termsimp(collect(int((a2+g2/sin(theta)^2)*p2*sin(theta),theta=0..Pi),lperts));

While this completes our formal analysis of the linear problem, these seven equations, must now be combined in some unobvious ways in order to facilitate the reduction of the source terms in the second order calculations.

> lsub1:={diff(k(r,t),t,t) = termsimp(collect(solve(expand((r-2*m)^2*leqn[5]+r^2*leqn[6]
-leqn[7]/r*(r-2*m)),diff(k(r,t),t,t)),lperts))};

In order to determine the second order PDEs ( qeqn_{1..7} ) we simply repeat the linear analysis, but now selecting only the epsilon^2 components of the Ricci tensor. The second order perturbation functions will, necessarily, appear in these PDEs in exactly the same manner as the linear perturbation functions appear in the leqn_i equations.

> aa1:=simplify(coeff(collect(grcomponent(R(dn,dn),[theta,theta]),epsilon),epsilon,2),
trig):

> gg1:=simplify(expand(1/sin(theta)^2*(coeff(collect(grcomponent(R(dn,dn),[phi,phi]),
epsilon),epsilon,2))),trig):

> zz1:=simplify(expand(gg1-aa1));

> zz2:=int(zz1*(2*cot(theta)*diff(p2,theta)+2*(2+1)*p2)*sin(theta),theta=0..Pi);

> qeqn1:= {H[2](r,t) = solve(zz2,H[2](r,t))};

> qeqn[1]:=qeqn1;

> grmap(R(dn,dn),subs,qeqn1,`x`);

Applying routine subs to R(dn,dn)

> aa1:=simplify(coeff(collect(grcomponent(R(dn,dn),[theta,theta]),epsilon),epsilon,2),
trig):

> gg1:=simplify(expand(1/sin(theta)^2*(coeff(collect(grcomponent(R(dn,dn),[phi,phi]),
epsilon),epsilon,2))),trig):

> qeqn[2]:=termsimp(collect(int((gg1+aa1)*p2*sin(theta),theta=0..Pi),qperts)):

> bb1:=simplify(coeff(collect(grcomponent(R(dn,dn),[theta,r]),epsilon),epsilon,2),trig):

> qeqn[3]:=termsimp(collect(int(bb1*diff(p2,theta)*sin(theta),theta=0..Pi),qperts)):

> cc1:=simplify(coeff(collect(grcomponent(R(dn,dn),[t,theta]),epsilon),epsilon,2),trig):

> qeqn[4]:=termsimp(collect(1/3*int(cc1*diff(p2,theta)*sin(theta),theta=0..Pi),qperts)):

> dd1:=simplify(coeff(collect(grcomponent(R(dn,dn),[r,t]),epsilon),epsilon,2),trig):

> qeqn[5]:=termsimp(collect(int(dd1*p2*sin(theta),theta=0..Pi),qperts)):

> ee1:=simplify(coeff(collect(grcomponent(R(dn,dn),[r,r]),epsilon),epsilon,2),trig):

> qeqn[6]:=termsimp(collect(expand(int(ee1*p2*sin(theta),theta=0..Pi)),qperts)):

> qsub1:={diff(H[0](r,t),t)=solve(qeqn[4],diff(H[0](r,t),t))};

> qeqn[1];

> undiff(qeqn[4]);

The derivation of the first Zerilli equation begins with the elimination of Hdot_0 using qeqn_4 , in the R^2_{rt} component of the Ricci tensor. We then make the previously defined Zerilli substitutions, and make various substitutions involving the leqn_i equations. (The actual substitutions are somewhat arbitrary. We choose to eliminate certain perturbation functions so as to reproduce the results presented in GL). The resulting expression is called etaeqn in our Maple worksheet. To simplify this rather awkward expression, we employ a new routine, ``mcollect()'' , that collects specified terms that are multiplied together. This routine could, for example, collect together terms that have a factor of, ab, but not, say, a^2, or b^2. A more complete description of this routine appears in Appendix A. (Readers who are familiar with the computer algebra environment will, no doubt, already see the utility of such a procedure.)

> eta1:=qeqn[5]:

> eta2:=termsimp(collect(subs(qsub1,eta1),qperts)):

> eta3:=termsimp(collect(simplify(subs(H[1](r,t)=newH(r,t),diff(K(r,t),t)=Kdot(r,t),
diff(K(r,t),r,t)=diff(Kdot(r,t),r),eta2)),zfuncs)):

> eta4:=eta(r,t)-termsimp(collect(solve(eta3,eta(r,t)),zfuncs)):

> eta5:=termsimp(collect(subs({diff(h[0](r,t),r,r)=solve(leqn[5],diff(h[0](r,t),r,r))},eta4),
zfuncs)):

> eta6:=termsimp(collect(subs({diff(h[0](r,t),r)=solve(leqn[2],diff(h[0](r,t),r)),
diff(k(r,t),r,t)=solve(leqn[4],diff(k(r,t),r,t))},eta5),zfuncs)):

> ug:={diff(k(r,t),r,r)=termsimp(collect(solve(subs(diff(h[0](r,t),r)=
solve(leqn[2],diff(h[0](r,t),r)),subs(lsub1,solve(leqn[6],diff(h[0](r,t),r,r))-
solve(leqn[5],diff(h[0](r,t),r,r)))),diff(k(r,t),r,r)),lperts))}:

> eta7:=termsimp(collect(subs(ug,eta6),zfuncs)):

> etaeqn:=eta7:

In order to produce a simplified final expression, we extract the source terms from etaeqn (by setting the second order perturbation functions to zero) and then simplify the result by repeated application of ``mcollect()'' .

> esource1:=mcollect(expand(eval(subs(eta(r,t)=0,chi(r,t)=0,etaeqn))), {h[1](r,t),h[0](r,t)}):

> esource2:=mcollect(esource1,{h[1](r,t),diff(h[1](r,t),t)}):

> esource3:=mcollect(esource2,{h[1](r,t),diff(k(r,t),r)}):

> esource4:=mcollect(esource3,{h[0](r,t),diff(k(r,t),t)}):

> esource5:=mcollect(esource4,{diff(h[1](r,t),t),diff(k(r,t),t)}):

> esource6:=mcollect(esource5,{k(r,t),h[1](r,t)}):

> esource7:=mcollect(esource6,{k(r,t),diff(k(r,t),t)}):

> esource8:=mcollect(esource7,{diff(k(r,t),r),diff(h[0](r,t),t)}):

> esource9:=mcollect(esource8,{k(r,t),diff(h[0](r,t),t)}):

> esource10:=mcollect(esource9,{h[0](r,t),diff(h[0](r,t),t)}):

> etasource:=kfactor(esource10,(r-2*m)/7/(2*r+3*m)):

It is now a simple matter to put these results together to arrive at the first Zerilli equation for the second order calculation:

> etafinal:={eta(r,t) = termsimp(collect(solve(eval(subs(h[0](r,t)=0,h[1](r,t)=0,k(r,t)=0,
etaeqn)),eta(r,t)),zfuncs))-etasource};

This agrees with the result found in [1].

The development of the second order Zerilli ``wave equation'' follows a similar procedure to that of the etafinal result. We first take the time derivative of qeqn_2 , eliminate Hdot_0 and Hdot^'_0 using qeqn_4 , make the Zerilli substitutions, and substitute for eta from our first Zerilli result. We then repeat this process for qenq_3 and add the resulting equations such that the chi^{..'} term vanishes. The result of these operations, which is called chieqn in our worksheet, is essentially the raw form of the result we are seeking. This expression is, however, rather large.

> atmp1:=diff(qeqn[2],t):

> atmp2:=termsimp(collect(expand(subs(qsub1,diff(qsub1,r),atmp1)),{diff(K(r,t),t,t,t)})):

> atmp3:=termsimp(collect(expand(eval(subs(H[1](r,t)=newH(r,t),
diff(K(r,t),t)=Kdot(r,t),diff(K(r,t),t,r)=diff(Kdot(r,t),r),diff(K(r,t),t,r,r)=diff(Kdot(r,t),r,r),
atmp2))),zfuncs)):

> atmp4:=termsimp(collect(expand(eval(subs(etafinal,atmp3))),zfuncs)):

> coefa:=simplify(coeff(atmp4,diff(chi(r,t),t,r,t),1));

> btmp1:=termsimp(collect(expand(subs(qsub1,diff(qsub1,r),diff(qeqn[3],t))),
{diff(H[1](r,t),t,t),diff(K(r,t),t,t,t)})):

> btmp2:=termsimp(collect(expand(eval(subs(H[1](r,t)=newH(r,t),diff(K(r,t),t)=Kdot(r,t),
diff(K(r,t),r,t)=diff(Kdot(r,t),r),btmp1))),zfuncs)):

> btmp3:=termsimp(collect(expand(eval(subs(etafinal,btmp2))),zfuncs)):

> coefb:=simplify(coeff(btmp3,diff(chi(r,t),r,t,t),1));

> ftmp1:=termsimp(collect(expand(atmp4/coefa-btmp3/coefb),{chi(r,t),diff(chi(r,t),r),
diff(chi(r,t),t),diff(chi(r,t),t,t),diff(chi(r,t),r,r),diff(chi(r,t),r,t),diff(chi(r,t),r,t,t),
diff(chi(r,t),r,r,r)})):

> ftmp2:=termsimp(collect(expand(eval(subs(k(r,t)=0,h[1](r,t)=0,h[0](r,t)=0,ftmp1))),
{chi(r,t),diff(chi(r,t),r),diff(chi(r,t),r,r),diff(chi(r,t),t,t),diff(chi(r,t),r,r,r)})):

> chi1:=termsimp(collect(diff(chi(r,t),t,t)-solve(ftmp1,diff(chi(r,t),t,t)),zfuncs)):

> chieqn:=ftmp1:

> nops(chieqn);

Although, at this point we can at least verify that it reproduces the correct linear Zerilli equation when the source terms are set to zero.

> termsimp(collect(expand(eval(subs(k(r,t)=0,h[1](r,t)=0,h[0](r,t)=0,
chieqn/coeff(chieqn,diff(chi(r,t),t,t),1)))),{chi(r,t),diff(chi(r,t),r),
diff(chi(r,t),r,r),diff(chi(r,t),t,t),diff(chi(r,t),r,r,r)}));

This is just the linear Zerilli wave equation with the quadrapole potential, as required.

In order to make contact with the results presented in GL, we now perform a ``gauge'' transformation and work with the new perturbation function zeta .

> zetasub:= {chi(r,t) = zeta(r,t)+2/7*(r^2/(2*r+3*m)*k(r,t)*diff(k(r,t),t)+k(r,t)^2)};

> zeta1:=eval(subs(zetasub,chi1)):

> nops(expand(zeta1));

This substitution has considerably reduced the number of terms in the expression.

> new2:=termsimp(collect(expand(leqn[4]-leqn[3]/2/m*3),lperts));

> leqn8:={diff(h[1](r,t),r) = termsimp(collect(solve(new2,diff(h[1](r,t),r)),lperts))};

> new3:=termsimp(collect(expand(solve(leqn[3],diff(h[1](r,t),r))-solve(diff(epxand(
leqn[4]/coeff(leqn[4],h[1](r,t),1)),r),diff(h[1](r,t),r))),lperts));

> leqn9:={h[1](r,t) = termsimp(collect(solve(new3,h[1](r,t)),lperts))};

> leqn10:={diff(h[1](r,t),t)=termsimp(collect(expand(solve(leqn[2],diff(h[1](r,t),t))),lperts))};

> aux1:=termsimp(collect(expand(solve(leqn[2],diff(h[1](r,t),t))-solve(diff(leqn[4],t),
diff(h[1](r,t),t))),lperts));

> aux2:=termsimp(collect(expand(diff((r-2*m)*r*leqn[4]/3,r))-leqn[3],lperts));

> aux3:=termsimp(collect(expand(solve(leqn[2],diff(h[1](r,t),t))-solve(leqn[7],
diff(h[1](r,t),t))),lperts));

> tmpaux:=termsimp(collect(expand(solve(leqn[5],diff(h[1](r,t),r,t))-solve(leqn[6],
diff(h[1](r,t),r,t))),lperts));

> tmpaux2:=termsimp(collect(expand(solve(tmpaux,diff(h[1](r,t),t))-solve(leqn[2],
diff(h[1](r,t),t))),lperts));

> tmpaux3:=termsimp(collect(expand(solve(tmpaux2,diff(k(r,t),t,t))-solve(leqn[7],
diff(k(r,t),t,t))),lperts));

> aux4:=termsimp(collect(expand(solve(tmpaux3,diff(h[1](r,t),t))-solve(leqn[2],
diff(h[1](r,t),t))),lperts));

> aux5:=termsimp(collect(solve(solve(diff(leqn[3],t),diff(h[1](r,t),r,t))-solve(leqn[5],
diff(h[1](r,t),r,t)),diff(h[1](r,t),t))-solve(leqn[2],diff(h[1](r,t),t)),lperts));

> aux6:=termsimp(collect(expand(solve(aux5,diff(k(r,t),t,t))-solve(aux3,diff(k(r,t),t,t))),
lperts));

Following GL, we now elimanate h_0 in favour of psi using:

> psisub:={h[0](r,t) = (psi(r,t) - r/3*k(r,t))*3*(2*r+3*m)/r/(r-2*m) +r*diff(k(r,t),r)};

This is the relation between the Zerilli function and the metric perturbation functions.

> aux7:=termsimp(collect(expand(subs(psisub,aux4)),lperts union {psi(r,t),
diff(psi(r,t),r)}));

> zsource1:=expand(simplify(subs(leqn8,zeta1))):

> nops(zsource1);

> zsource2:=expand(simplify(eval(subs(leqn10,zsource1)))):

> nops(zsource2);

> zsource3:=expand(simplify(subs(leqn9,zsource2))):

> nops(zsource3);

> zsource3a:=expand(simplify(eval(subs(diff(k(r,t),t,t,r)=solve(aux1,diff(k(r,t),t,t,r)),
zsource3)))):

> zsource4:=expand(simplify(eval(subs(diff(k(r,t),t,t)=solve(aux3,diff(k(r,t),t,t)),
zsource3a)))):

> zsource5:=expand(simplify(eval(subs(diff(k(r,t),r,r,t)=solve(aux2,diff(k(r,t),r,r,t)),
zsource4)))):

> zsource6:=expand(simplify(eval(subs(diff(h[0](r,t),t,t)=solve(aux6,diff(h[0](r,t),t,t)),
zsource5)))):

> nops(zsource6);

> zsource7:=expand(simplify(eval(subs(psisub,zsource6)))):

> nops(zsource7);

> zsource8:=expand(simplify(eval(subs(k(r,t)=solve(aux7,k(r,t)),zsource7)))):

> zetaeqn:=zsource8:

> s_ren := make_s():

> sprime:=make_s_mu():

After making this substitution, and a fairly extensive simplification of the source terms using the linear perturbation equations, we arrive at the unsimplified answer zetaeqn.

> nops(zetaeqn);

Here again we make use of the ``mcollect()'' routine to simplify our result.

> final1:=mcollect(expand(zetaeqn),{psi(r,t),diff(psi(r,t),r,r)}):

> final2:=mcollect(final1,{psi(r,t),diff(psi(r,t),t)}):

> final3:=mcollect(final2,{diff(psi(r,t),t),diff(psi(r,t),r,r,r)}):

> final4:=mcollect(final3,{diff(psi(r,t),r),diff(psi(r,t),t)}):

> final5:=mcollect(final4,{diff(psi(r,t),t),diff(psi(r,t),t,r)}):

> final6:=mcollect(final5,{psi(r,t),diff(psi(r,t),r)}):

> final7:=mcollect(final6,{diff(psi(r,t),r),diff(psi(r,t),t,r,r)}):

> final8:=mcollect(final7,{diff(psi(r,t),t),diff(psi(r,t),r,r)}):

> final9:=mcollect(final8,{psi(r,t),diff(psi(r,t),r,t)}):

> final10:=mcollect(final9,{diff(psi(r,t),r),diff(psi(r,t),r,t)}):

> final11:=mcollect(final10,{psi(r,t),diff(psi(r,t),t,r,r)}):

> final12:=mcollect(final11,{diff(psi(r,t),r),diff(psi(r,t),r,r)}):

> final13:=mcollect(final12,{psi(r,t)}):

> final14:=mcollect(final13,{diff(psi(r,t),t)}):

> final15:=mcollect(final14,{diff(psi(r,t),r)}):

> final16:=mcollect(final15,{diff(psi(r,t),r,r)}):

> final17:=mcollect(final16,{diff(psi(r,t),r,t)}):

> nops(final17);

We now make one final set of substitutions, mu=(r-2m) , and, lambda=(2r+3m) , and extract the source terms from the second order Zerilli wave equation:

> zetasource:=subs((r-2*m)=mu,(2*m-r)=-mu,(2*r+3*m)=lambda,
kfactor(eval(subs(zeta(r,t)=0,-final17)),12/7*(r-2*m)^3/(2*r+3*m)));

This is our final result.

> final18:=subs((r-2*m)=mu,(-r+2*m)=-mu,(2*r+3*m)=lambda,final17):

> final19:=kfactor(collect(final18,{zeta(r,t),diff(zeta(r,t),r),diff(zeta(r,t),r,r),
diff(zeta(r,t),t,t)}),12/7*mu^3/lambda):

> simplify(expand(final19-final18));

> termsimp(simplify(expand(subs(zeta(r,t)=0,lambda=(2*r+3*m),mu=(r-2*m),
final19)+s_ren)));

> termsimp(simplify(expand(subs(zeta(r,t)=0,lambda=(2*r+3*m),mu=(r-2*m),final19)
+sprime)));

>