Linear perturbation of the Schwarzschild metric

Preliminaries

Kernels for parallel evaluation

In the following we will represent tensor fields by arrays, i.e., nested lists. To reduce the time needed for simplification, we simplify their components in parallel. The following command launches several parallel kernels, corresponding to the number of available CPU cores.

LaunchKernels[];

Assumptions on coordinates

To aid with the simplification, we specify the ranges of the coordinates and parameters we use.

$Assumptions={r>0,-1≤χ≤1,0≤φ≤2Pi,m∈Integers,l∈Integers};

ParallelEvaluate[$Assumptions={r>0,-1≤χ≤1,0≤φ≤2Pi,m∈Integers,l∈Integers}];

Coordinates

The coordinates we use. Note that we use χ=cos ϑ instead of the azimuth angle ϑ.

co={t,r,χ,φ}

Simplification to linear order

Since we are working with linear perturbations, which we define with respect to a perturbation parameter ϵ, we discard all terms of higher than linear order, and then simplify.

linsimp[x_]:=FullSimplify[Normal[Series[x,{ϵ,0,1}]]]

For convenience, we also define a function with simplifies all tensor components in parallel.

tensimp[x_]:=ParallelMap[linsimp,x,{TensorRank[x]}]

Formatting

We will use a number of functions to express the metric perturbations. In order to shorten the expressions, we define a number of shorthand notations for these functions and their derivatives.

Format[Y[l_,m_,χ_,φ_]]:=Subscript["Y",StringJoin@@(ToString/@{l,m})]

Format[Derivative[0,0,c_,0][Y][l_,m_,χ_,φ_]]:=Subsuperscript["Y",StringJoin@@(ToString/@{l,m}),StringRepeat["'",c]]

Format[sp1[t_,r_]]:=Subscript["S","tt"]

Format[sp2[t_,r_]]:=Subscript["S","tr"]

Format[sp3[t_,r_]]:=Subscript["S","rr"]

Format[vp1[t_,r_]]:=Subscript["V","t"]

Format[vp2[t_,r_]]:=Subscript["V","r"]

Format[va1[t_,r_]]:=Subscript["A","t"]

Format[va2[t_,r_]]:=Subscript["A","r"]

Format[tt[t_,r_]]:="T"

Format[tp[t_,r_]]:="P"

Format[ta[t_,r_]]:="E"

Format[Derivative[ct_,cr_][sp1][t_,r_]]:=Subscript["S",StringJoin["tt,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][sp2][t_,r_]]:=Subscript["S",StringJoin["tr,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][sp3][t_,r_]]:=Subscript["S",StringJoin["rr,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][vp1][t_,r_]]:=Subscript["V",StringJoin["t,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][vp2][t_,r_]]:=Subscript["V",StringJoin["r,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][va1][t_,r_]]:=Subscript["A",StringJoin["t,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][va2][t_,r_]]:=Subscript["A",StringJoin["r,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][tt][t_,r_]]:=Subscript["T",StringJoin[",",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][tp][t_,r_]]:=Subscript["P",StringJoin[",",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][ta][t_,r_]]:=Subscript["E",StringJoin[",",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[ν[r_]]:="ν"

Format[λ[r_]]:="λ"

Format[Derivative[c_][ν][r_]]:=Superscript["ν",StringRepeat["'",c]]

Format[Derivative[c_][λ][r_]]:=Superscript["λ",StringRepeat["'",c]]

Format[ξt[t_,r_]]:=Subscript["ξ","t"]

Format[ξr[t_,r_]]:=Subscript["ξ","r"]

Format[ξp[t_,r_]]:="Π"

Format[ξa[t_,r_]]:="Ξ"

Format[Derivative[ct_,cr_][ξt][t_,r_]]:=Subscript["ξ",StringJoin["t,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][ξr][t_,r_]]:=Subscript["ξ",StringJoin["r,",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][ξp][t_,r_]]:=Subscript["Π",StringJoin[",",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[Derivative[ct_,cr_][ξa][t_,r_]]:=Subscript["Ξ",StringJoin[",",StringRepeat["t",ct],StringRepeat["r",cr]]]

Format[fsp1[r_]]:=Subscript["(OverHat[S])","tt"]

Format[fsp2[r_]]:=Subscript["(OverHat[S])","tr"]

Format[fsp3[r_]]:=Subscript["(OverHat[S])","rr"]

Format[fva1[r_]]:=Subscript["(OverHat[A])","t"]

Format[fva2[r_]]:=Subscript["(OverHat[A])","r"]

Format[ftt[r_]]:="(OverHat[T])"

Format[Derivative[c_][fsp1][r_]]:=Subsuperscript["(OverHat[S])","tt",StringRepeat["'",c]]

Format[Derivative[c_][fsp2][r_]]:=Subsuperscript["(OverHat[S])","tr",StringRepeat["'",c]]

Format[Derivative[c_][fsp3][r_]]:=Subsuperscript["(OverHat[S])","rr",StringRepeat["'",c]]

Format[Derivative[c_][fva1][r_]]:=Subsuperscript["(OverHat[A])","t",StringRepeat["'",c]]

Format[Derivative[c_][fva2][r_]]:=Subsuperscript["(OverHat[A])","r",StringRepeat["'",c]]

Format[Derivative[c_][ftt][r_]]:=Superscript["(OverHat[T])",StringRepeat["'",c]]

Format[fax[r_]]:="𝒜"

Format[Derivative[c_][fax][r_]]:=Superscript["𝒜",StringRepeat["'",c]]

Two-dimensional geometry on the sphere

Metric

We define the metric (Subscript[γ, ab]) in the angular coordinates (χ,φ) we use.

MatrixForm[Met2DD=DiagonalMatrix[{1/(1-χ^2),1-χ^2}]]

Inverse metric

The inverse metric is immediate.

MatrixForm[Met2UU=Inverse[Met2DD]]

Antisymmetric tensor

We define the totally antisymmetric tensor on the sphere with a negative sign, since the coordinate frame in the (χ,φ) has the opposite orientation compared to the standard (ϑ,φ) coordinates, and we use the latter orientation.

MatrixForm[Eps2DD=-LeviCivitaTensor[2]Sqrt[Det[Met2DD]]]

Christoffel symbols

Finally, the Levi-Civita connection of the metric (Subscript[γ, ab]).

TableForm[Gam2UDD=Table[Sum[Met2UU[[a,d]](D[Met2DD[[d,c]],co[[b+2]]]+D[Met2DD[[b,d]],co[[c+2]]]-D[Met2DD[[b,c]],co[[d+2]]])/2,{d,2}],{a,2},{b,2},{c,2}]]

Tensor harmonics

Derivatives of spherical harmonics

Here we implement the rules for derivatives of the spherical harmonics (Subscript[Y, lm]) by hand, since we use rather non-standard coordinates. These are defined by the Laplace equation and the usual angular dependence.

yru={Derivative[0,0,2,0][Y][l_,m_,χ_,φ_]:>(2Derivative[0,0,1,0][Y][l,m,χ,φ]χ-(l(l+1)-m^2/(1-χ^2))Y[l,m,χ,φ])/(1-χ^2),Derivative[0,0,a_,b_/;b>0][Y][l_,m_,χ_,φ_]:>Derivative[0,0,a,0][Y][l,m,χ,φ](I m)^b};

This is their effect.

D[Y[l,m,χ,φ],{χ,2}]/.yru

D[Y[l,m,χ,φ],φ]/.yru

We can also define them explicitly.

YY[l_,m_,χ_,φ_]:=Sqrt[(2l+1)(l-m)!/(l+m)!/4/Pi]LegendreP[l,m,χ]Exp[I m φ]

These are the first few harmonics.

TableForm[Table[YY[l,m,χ,φ],{l,0,2},{m,-l,l}]]

Scalar harmonic

The scalar harmonics are simply the spherical harmonics themselves. They have polar parity.

YS[l_,m_,χ_,φ_]:=Y[l,m,χ,φ]

The first non-vanishing scalar is the monopole l = 0.

YS[0,0,χ,φ]/.Y->YY

Vector harmonics

For the vector harmonics, we have two contributions. The first one has polar parity.

YVP[l_,m_,χ_,φ_]:=D[YS[l,m,χ,φ],#]&/@{χ,φ}//.yru

YVP[l,m,χ,φ]

The second part has axial parity, which can be seen from the appearance of the Levi-Civita tensor.

YVA[l_,m_,χ_,φ_]:=Eps2DD.Met2UU.YVP[l,m,χ,φ]

YVA[l,m,χ,φ]

For the vectors, the monopole vanishes, and the first non-vanishing term is the dipole l = 1.

TableForm[Map[MatrixForm,Table[FullSimplify[YVP[l,m,χ,φ]/.Y->YY],{l,0,1},{m,-l,l}],{2}]]

TableForm[Map[MatrixForm,Table[FullSimplify[YVA[l,m,χ,φ]/.Y->YY],{l,0,1},{m,-l,l}],{2}]]

Tensor harmonics

Finally, we come to the tensor harmonics. The first one is just the diagonal trace component. It has polar parity.

YTT[l_,m_,χ_,φ_]:=YS[l,m,χ,φ]Met2DD

MatrixForm[YTT[l,m,χ,φ]]

Also the second covariant derivative gives a term with polar parity.

YTP[l_,m_,χ_,φ_]:=Module[{a,b,c},Table[D[YS[l,m,χ,φ],co[[a+2]],co[[b+2]]]-Sum[D[YS[l,m,χ,φ],co[[c+2]]]Gam2UDD[[c,a,b]],{c,2}],{a,2},{b,2}]//.yru]

MatrixForm[FullSimplify[YTP[l,m,χ,φ]]]

The third and last term has axial parity.

YTA[l_,m_,χ_,φ_]:=(#+Transpose[#])/2&[Eps2DD.Met2UU.YTP[l,m,χ,φ]]

MatrixForm[FullSimplify[YTA[l,m,χ,φ]]]

Note that the axial tensor has vanishing trace by definition:

Tr[Met2UU.YTA[l,m,χ,φ]]

The trace of the polar tensor, however, is given by the Laplacian:

FullSimplify[Tr[Met2UU.YTP[l,m,χ,φ]]]

Often one therefore also uses a trace-free polar tensor instead:

FullSimplify[Tr[Met2UU.(YTP[l,m,χ,φ]+l(l+1)YTT[l,m,χ,φ]/2)]]

For the two trace-free tensor harmonics, the lowest order is the quadrupole l = 2.

TableForm[Map[MatrixForm,Table[FullSimplify[YTP[l,m,χ,φ]+l(l+1)YTT[l,m,χ,φ]/2/.Y->YY],{l,0,2},{m,-l,l}],{2}]]

TableForm[Map[MatrixForm,Table[FullSimplify[YTA[l,m,χ,φ]/.Y->YY],{l,0,2},{m,-l,l}],{2}]]

Calculate Laplacian

It is instructive to calculate the Laplacian of the tensor harmonics, which is one possibility to define them (along with the action of the angular momentum operators).

FullSimplify[Sum[Met2UU[[a,b]](D[YS[l,m,χ,φ],co[[a+2]],co[[b+2]]]-Sum[D[YS[l,m,χ,φ],co[[c+2]]]Gam2UDD[[c,a,b]],{c,2}]),{a,2},{b,2}]/YS[l,m,χ,φ]//.yru]

Spacetime geometry

Background metric

We proceed with the spacetime geometry. This is the unperturbed background metric. Here we assume a spherically symmetric, static metric, but we leave the (Subscript[g, tt]) and (Subscript[g, rr]) components as free functions of the radial coordinate for now.

MatrixForm[Met4DD0=DiagonalMatrix[{-Exp[2ν[r]],Exp[2λ[r]],r^2/(1-χ^2),r^2(1-χ^2)}]]

Metric perturbation

We then define the metric perturbation.

MatrixForm[Met4DD1=tensimp[ArrayFlatten[{{{{sp1[t,r]YS[l,m,χ,φ]}},{{sp2[t,r]YS[l,m,χ,φ]}},{vp1[t,r]YVP[l,m,χ,φ]+va1[t,r]YVA[l,m,χ,φ]}},{{{sp2[t,r]YS[l,m,χ,φ]}},{{sp3[t,r]YS[l,m,χ,φ]}},{vp2[t,r]YVP[l,m,χ,φ]+va2[t,r]YVA[l,m,χ,φ]}},{Transpose[{vp1[t,r]YVP[l,m,χ,φ]+va1[t,r]YVA[l,m,χ,φ]}],Transpose[{vp2[t,r]YVP[l,m,χ,φ]+va2[t,r]YVA[l,m,χ,φ]}],tt[t,r]YTT[l,m,χ,φ]+tp[t,r]YTP[l,m,χ,φ]+ta[t,r]YTA[l,m,χ,φ]}}]//.yru]]

Full metric

The full metric is given as the sum of the background metric and the perturbation. Note that we introduce the perturbation parameter ϵ here, which will be used to select only linear perturbations by performing a Taylor expansion.

Met4DD=Met4DD0+ϵ Met4DD1;

Inverse metric

This is the inverse metric. Note that only terms up to linear order in ϵ appear, since we have discarded all higher order terms in the simplification function.

MatrixForm[Met4UU=tensimp[Inverse[Met4DD]]]

Christoffel symbols

Next, we calculate the Christoffel symbols.

Gam4UDD=tensimp[Table[Sum[Met4UU[[a,d]](D[Met4DD[[d,c]],co[[b]]]+D[Met4DD[[b,d]],co[[c]]]-D[Met4DD[[b,c]],co[[d]]])/2,{d,4}],{a,4},{b,4},{c,4}]//.yru];

It is helpful to note also the background of the Christoffel symbols. It can be used to calculate the covariant derivative of the metric perturbation, which then yields the perturbation of the Riemann tensor.

TableForm[FullSimplify[Gam4UDD/.ϵ->0]]

Riemann tensor

This gives us the perturbed Riemann tensor.

Riem4UDDD=tensimp[Table[D[Gam4UDD[[a,b,d]],co[[c]]]-D[Gam4UDD[[a,b,c]],co[[d]]]+Sum[Gam4UDD[[a,e,c]]Gam4UDD[[e,b,d]]-Gam4UDD[[a,e,d]]Gam4UDD[[e,b,c]],{e,4}],{a,4},{b,4},{c,4},{d,4}]//.yru];

Also the background value of the Riemann tensor is easily obtained.

TableForm[FullSimplify[Riem4UDDD/.ϵ->0]]

Ricci tensor

Contracting two indices, we get the Ricci tensor.

Ric4DD=tensimp[TensorContract[Riem4UDDD,{{1,3}}]];

This is the background of the Ricci tensor.

TableForm[Ric4DD0=tensimp[Ric4DD/.ϵ->0]]

For later use, we also isolate the perturbation from the background.

Ric4DD1=tensimp[D[Ric4DD,ϵ]];

Check spherical symmetry

Killing vector fields on the sphere

We can check the expansion of the perturbations into tensor spherical harmonics by calculating their behavior under rotations. For this purpose, we need the rotation generating vector fields on the sphere, in the relevant coordinates.

kvf2=tensimp[{{Sin[φ],Cos[φ]/Tan[ϑ]},{-Cos[φ],Sin[φ]/Tan[ϑ]},{0,-1}}.DiagonalMatrix[{-Sin[ϑ],1}]/.{ϑ->ArcCos[χ]}]

We see that they indeed satisfy the rotation algebra.

tensimp[Table[Sum[kvf2[[i,b]]D[kvf2[[j,a]],co[[b+2]]]-kvf2[[j,b]]D[kvf2[[i,a]],co[[b+2]]],{b,2}]-Sum[Signature[{i,j,k}]kvf2[[k,a]],{k,3}],{i,3},{j,3},{a,2}]]

Killing vector fields on spacetime

We lift the vector fields to the full spacetime manifold by noticing that they do not modify the time and radial coordinates.

kvf4=PadLeft[#,4]&/@kvf2

Symmetry of two-dimensional geometry

By construction, the two-dimensional metric and its Levi-Civita tensor on the sphere are invariant under rotations.

tensimp[Table[Sum[kvf2[[i,c]]D[#[[a,b]],co[[c+2]]]+#[[c,b]]D[kvf2[[i,c]],co[[a+2]]]+#[[a,c]]D[kvf2[[i,c]],co[[b+2]]],{c,2}],{i,3},{a,2},{b,2}]&/@{Met2DD,Eps2DD}]

Symmetry of four-dimensional background geometry

Also the spacetime metric is invariant under rotations by construction.

tensimp[Table[Sum[kvf4[[i,c]]D[Met4DD0[[a,b]],co[[c]]]+Met4DD0[[c,b]]D[kvf4[[i,c]],co[[a]]]+Met4DD0[[a,c]]D[kvf4[[i,c]],co[[b]]],{c,4}],{i,3},{a,4},{b,4}]]

Transformation of perturbation

Finally, we calculate how the metric perturbations transform under rotations.

rotpert=tensimp[Table[Sum[kvf4[[i,c]]D[Met4DD1[[a,b]],co[[c]]]+Met4DD1[[c,b]]D[kvf4[[i,c]],co[[a]]]+Met4DD1[[a,c]]D[kvf4[[i,c]],co[[b]]],{c,4}],{i,3},{a,4},{b,4}]//.yru];

We see that rotation around the z axis simply gives us the magnetic quantum number m.

TableForm[tensimp[rotpert[[3]]+I m Met4DD1]]

Rotation around the other axes corresponds to the ladder operators. Hence, they allow us to change the value of m.

TableForm[tensimp[rotpert[[1]]+rotpert[[2]]I-(Met4DD1/.{Y->(Exp[I#4]I(Derivative[0,0,1,0][Y][#1,#2-1,#3,#4]Sqrt[1-#3^2]+Y[#1,#2-1,#3,#4](#2-1)#3/Sqrt[1-#3^2])&)}//.yru/.{m->m+1})]]

TableForm[tensimp[rotpert[[1]]-rotpert[[2]]I-(Met4DD1/.{Y->(Exp[-I#4]I(-Derivative[0,0,1,0][Y][#1,#2+1,#3,#4]Sqrt[1-#3^2]+Y[#1,#2+1,#3,#4](#2+1)#3/Sqrt[1-#3^2])&)}//.yru/.{m->m-1})]]

Regge-Wheeler gauge

Gauge-transforming vector field

We can simplify the following calculations by assuming a particular gauge. To arrive at this gauge, we first consider a general vector field, at the same order as the metric perturbations, and expand it in spherical harmonics.

TableForm[rgvf=Inverse[Met4DD0].Join[{ξt[t,r],ξr[t,r]}YS[l,m,χ,φ],(ξp[t,r]YVP[l,m,χ,φ]+ξa[t,r]YVA[l,m,χ,φ])/r^2]]

Change of the metric perturbation

The change of the linear metric perturbation is given by the Lie derivative of the background metric with respect to the transforming vector field.

TableForm[pertchg=tensimp[Table[Sum[rgvf[[c]]D[Met4DD0[[a,b]],co[[c]]]+Met4DD0[[c,b]]D[rgvf[[c]],co[[a]]]+Met4DD0[[a,c]]D[rgvf[[c]],co[[b]]],{c,4}],{a,4},{b,4}]//.yru]]

Transformed perturbation

To study how the different components of the metric change, we now calculate their sum - this is the perturbation of the metric after performing the gauge transformation.

gtmet=tensimp[pertchg+Met4DD1];

Decomposition of components

In order to choose a particularly convenient gauge for the further calculations, it is helpful to study how the different components of the perturbation transform. For this purpose, it is helpful to first gather the time and radial components of the gauge transforming covector field into a two-dimensional object.

xitr={ξt[t,r],ξr[t,r]}

The three scalar components transform as a symmetric two-tensor in the time and radial coordinates. Hence, we get their transformation from a symmetrized covariant derivative of the two-dimensional covector field.

gtsca=Union[Flatten[MapThread[Rule,{{{sp1[t,r],sp2[t,r]},{sp2[t,r],sp3[t,r]}},-Table[D[xitr[[i]],co[[j]]]+D[xitr[[j]],co[[i]]]-2Sum[Gam4UDD[[k,i,j]]xitr[[k]],{k,2}],{i,2},{j,2}]/.ϵ->0},2]]]

The axial vector components transform as the gradient of the axial component of the gauge vector field, properly renormalized with the radius.

gtaxi=MapThread[Rule[#1,-D[ξa[t,r]/r^4,#2]r^2]&,{{va1[t,r],va2[t,r]},{t,r}},1]

Similarly, the polar vector components transform as the gradient of the polar component of the gauge vector field, but with an additional contribution from the time and radial part.

gtvec=MapThread[Rule[#1,-#2-D[ξp[t,r]/r^4,#3]r^2]&,{{vp1[t,r],vp2[t,r]},xitr,{t,r}},1]

Finally, the tensor components have a simple transformation behavior.

gtten={tp[t,r]->-2ξp[t,r]/r^2,ta[t,r]->-2ξa[t,r]/r^2,tt[t,r]->-2Exp[-2λ[r]]ξr[t,r]r}

We check that these indeed give the full transformation.

TableForm[tensimp[gtmet/.gtsca/.gtaxi/.gtvec/.gtten]]

Gauge choice

From the gauge transformations above, we see that by an appropriate choice of the gauge vector field we can eliminate the two off-diagonal tensor components of the perturbation, as well as the polar vector parts. This is achieved by choosing the following vector field:

rwgauge={ξp->(-#2^2tp[#1,#2]/2&),ξa->(-#2^2ta[#1,#2]/2&),ξt->(Derivative[1,0][tp][#1,#2]/2-vp1[#1,#2]&),ξr->(Derivative[0,1][tp][#1,#2]/2-tp[#1,#2]/#2-vp2[#1,#2]&)}

We see that indeed now the tensor part becomes diagonal, and the vector part retains only axial components:

TableForm[rwmet=tensimp[gtmet/.rwgauge]]

Hence, we can assume that the metric perturbation is already given in this gauge, which is called the Regge-Wheeler gauge.

rwru={ta->(0&),tp->(0&),vp1->(0&),vp2->(0&)};

One finds that imposing the Regge-Wheeler gauge on either the transformed or untransformed metric perturbation gives the same result, since it means that the gauge transformation becomes trivial if one already starts from the Regge-Wheeler gauge.

TableForm[tensimp[Met4DD1/.rwru]]

TableForm[tensimp[rwmet/.rwru]]

Ricci tensor perturbation

For later use, we now calculate the Ricci tensor perturbation in the Regge-Wheeler gauge.

Ric4DD1RW=tensimp[Ric4DD1/.rwru];

Perturbation of the Schwarzschild metric

Schwarzschild metric

In the following, we will assume that the background solution is given by the Schwarzschild metric.

sbhru={ν->(Log[1-R/#]/2&),λ->(-Log[1-R/#]/2&)};

MatrixForm[Met4DD0/.sbhru]

Background solution

We check that this is indeed a vacuum solution of Einstein’s equations. Of course, this is known - here we just double check.

tensimp[Ric4DD0/.sbhru]

Ricci tensor perturbation

We then evaluate the perturbation of the Ricci tensor for the Schwarzschild metric.

Ric4DD1RWS=tensimp[Ric4DD1RW/.sbhru];

Fourier transformation

In the next step, we perform a Fourier transformation with respect to the time coordinate. Here we consider only a single Fourier mode with frequency ω. We thus separate the perturbations into time and radial dependence.

fourier={tt->(Exp[I ω#1]ftt[#2]&),sp1->(Exp[I ω#1]fsp1[#2]&),sp2->(Exp[I ω#1]fsp2[#2]&),sp3->(Exp[I ω#1]fsp3[#2]&),va1->(Exp[I ω#1]fva1[#2]&),va2->(Exp[I ω#1]fva2[#2]&)};

In the resulting formula for the Ricci tensor perturbation, we find factors iω in place of the time derivatives.

Ric4DD1RWSF=tensimp[Exp[-I ω t]Ric4DD1RWS/.fourier];

Case m=0

We finally restrict ourselves to the case m=0, since the spherical symmetry implies that the solution will not depend on m.

Ric4DD1RWSF0=tensimp[Ric4DD1RWSF/.{m->0}];

Decomposing the axial equations

Non-trivial axial equations

Looking at the perturbation of the Ricci tensor, one can see that the axial and polar perturbations decouple. Here we focus on the axial equations. We see that there are three equations, which after canceling some common factors are given as follows.

eq03=FullSimplify[2r^3Ric4DD1RWSF0[[1,4]]/(-1+χ^2)/D[Y[l,0,χ,φ],χ]]

eq13=FullSimplify[2r^2(r-R)Ric4DD1RWSF0[[2,4]]/(-1+χ^2)/D[Y[l,0,χ,φ],χ]]

eq23=FullSimplify[2r^2(r-R)Ric4DD1RWSF0[[3,4]]/(l(l+1)Y[l,0,χ,φ]-2D[Y[l,0,χ,φ],χ]χ)]

Relation between equations

Note first that the three equations above are not independent of each other. It turns out that the first equation is a linear combination of the second, the third and the derivative of the second equations. To see this, we gather these four equations into a list.

TableForm[eqnax={eq03,eq13,eq23,D[eq13,r]}]

The independent functions which appear in these equations are the following axial perturbation components and their derivatives.

terms={fva1[r],fva1'[r],fva1''[r],fva2[r],fva2'[r]}

Now we have four equations and five terms which appear in them. The coefficients of these five terms are given as follows.

TableForm[coeff=Outer[Coefficient,eqnax,terms]]

If the equations are independent, then the matrix above has maximal rank, which is 4. However, we see that this is not the case.

MatrixRank[coeff]

It follows that the transposed matrix has a non-trivial, one-dimensional kernel, which we can find as follows (and renormalize for simplicity).

ns=NullSpace[Transpose[coeff]][[1]]r(r-R)

Finally, we check that the corresponding linear combination of the equations indeed vanishes.

FullSimplify[eqnax.ns]

Solve equation for time component

We are now left with two equations for two unknown functions. Note that the last equation does not contain derivatives acting on the time component of the perturbation. Hence, we may solve it algebraically.

sol23={fva1->(Evaluate[FullSimplify[Solve[eq23==0,fva1[r]][[1]]][[1,2]]/.{r->#1}]&)}

Remaining equation

To obtain an equation for the radial component, we insert the solution for the time component into the only remaining equation.

eqax=FullSimplify[eq13/.sol23]

Simplification using tortoise coordinate

To simplify the last equation, we change the radial coordinate. The coordinate we choose here is called the tortoise coordinate. We use the following transformation.

tortoise={fva2->(fax[#+Log[#/R-1]R]#/(1-R/#)&)};

This is the final equation we obtain. Note that the appearing derivative is with respect to the tortoise coordinate, not with respect to the original radial coordinate.

eqaxt=FullSimplify[(r-R)eqax/r^5/.tortoise]