TextCell[
"Even Parity Perturbations of the Self-Similar LTB Spacetime", "Title" ]

TextCell["In this notebook we study linear even parity perturbations of the self-similar LTB spacetime. Section 5.2 of the thesis discusses the initial system of equations and the procedure we follow in reducing them. The details of this procedure are provided here and Section 5.2 of the thesis should be read in conjunction with this notebook. 

In this notebook, we first input the Gerlach-Sengupta perturbations and the linearised Einstein equations. We introduce a first order reduction and identify the initial set of perturbation equations. This set consists of seven equations for six variables. We identify the Einstein constraint, and the trivial constraint which arises from the first order reduction. We then use these two constraints to eliminate two variables, thus reducing the initial system to a four dimensional system (that is, four equations in four variables). We do this in two stages, first reducing to a five dimensional, and then to a four dimensional system. In so doing, we identify the five dimensional symmetric hyperbolic system which is used in Chapters 5 and 6 of the thesis. 

The notebook is divided into sections as follows: 
Section 1: The Background Spacetime. In this section, we input the background metric and calculate the connection components and the Ricci tensor.
Section 2: The Perturbation Variables and Equations of Motion. Here we input the Gerlach-Sengupta equations and the perturbation terms for the even parity perturbations. 
Section 3: The Initial System. In this section we take the initial equations and rewrite them as a vector equation with matrix coefficients. We identify two constraints, namely the Einstein constraint and the trivial constraint.  
Section 4: Constraint Propagation. In this section, we show explicitly that the Einstein constraint is propagated by the equations of motion. 
Section 5: Dimensional Reduction. In this section, we reduce the system first to a six dimensional, and then to a five dimensional system. 
Section 6: The Four Dimensional System. We eliminate the final constraint and reduce to the four dimensional system. We note that this system is not symmetric hyperbolic. 
   ", "Text", CellFrame->True, Background->GrayLevel[0.9]]

Clear["Global`*"]

x={z,p,theta,phi}; 

(* We input the background metric, with the rr component written in terms of S only. See Chapter 2, equation 2.9 *)
gdd=Exp[2*p]*{{-1,-z,0,0},{-z,(-z^2+(S[z]-z*D[S[z], z])^2),0,0},{0,0,(S[z])^2,0},{0,0,0,(S[z])^2*Sin[theta]^2}};

(* The metric with indices upstairs: *)
guu=Simplify[Inverse[gdd]];

(* The connection components: *)
For[i=1,i<5,i++,For[j=i, j<5, j++,For[k=1, k<5, k++, 
			chrisddd[j,i,k]=chrisddd[i,j,k]= (D[ gdd[[j,k]],x[[i]] ]+D[ gdd[[i,k]],x[[j]] ]-D[ gdd[[i,j]],x[[k]]])/2 ]]];

For[i=1,i<5,i++, 
	For[j=i, j<5, j++,
				For[k=1,k<5,k++, chrisddu[j,i,k]=chrisddu[i,j,k]=Sum[guu[[k,l]]*chrisddd[i,j,l],{l,1,4}] ]]];

(* The Ricci tensor: *)
For[i=1,i<5,i++,
	For[j=i,j<5,j++,
		ric1dd[j,i]=ric1dd[i,j]=Sum[ (D[chrisddu[i,j,n],x[[n]]]- D[chrisddu[j,n,n],x[[i]]]), {n,1,4}] ]];
For[i=1,i<5,i++,
	For[j=i,j<5,j++,
		 ric2dd[j,i]=ric2dd[i,j]=Sum[ Sum[(chrisddu[k,n,n]*chrisddu[i,j,k]- chrisddu[i,k,n]*chrisddu[n,j,k]) , {k,1,4}], {n,1,4}] ]];
For[i=1,i<5,i++,
	For[j=i,j<5,j++,  ric4dd[j,i]=ric4dd[i,j]=Simplify[ric1dd[i,j]+ric2dd[i,j]]]];

For[i=1,i<5,i++,
	For[j=i,j<5,j++,  ric4uu[j,i]=ric4uu[i,j]=Sum[Sum[guu[[i,k]]*guu[[j,l]]*ric4dd[k,l],{k,1,4}],{l,1,4}]]];
ricci=Sum[guu[[i,j]]ric4dd[i,j],{i,1,4},{j,1,4}];

sform[z_]=(1+a*z)^(2/3);
qform[z_]=2*a^2*(1/(4*9*Pi*S[z]^2*(S[z]-z*D[S[z], z])));
S[z_]=sform[z];

Simplify[ricci]
Simplify[ricci-8*Pi*qform[z]*Exp[-2p]]
Clear[ S]

Clear[x,gdd,guu,chrisddd,chrisddu,ric1dd,ric2dd,ricdd];

(* The coordinates on M^2: *)
x={z,p};

(* The metric on M^2.  *)
gdd=Exp[2*p]*{{-1,-z},{-z,(-z^2+(S[z]-z*D[S[z], z])^2)}};

guu=Simplify[Inverse[gdd]];

(* The connection components: *)
For[i=1,i<3,i++, 
	For[j=i, j<3, j++,For[k=1, k<3, k++, 
			chrisddd[j,i,k]=chrisddd[i,j,k]= (D[ gdd[[j,k]],x[[i]] ]+D[ gdd[[i,k]],x[[j]] ]- D[ gdd[[i,j]],x[[k]]])/2 ]]];
For[i=1,i<3,i++, 
	For[j=i, j<3, j++,
				For[k=1,k<3,k++,
				chrisddu[j,i,k]=chrisddu[i,j,k]= Sum[guu[[k,l]]*chrisddd[i,j,l],{l,1,2}] ]]];

(* The Ricci tensor: *)
For[i=1,i<3,i++,
	For[j=i,j<3,j++,
		ric1dd[j,i]=ric1dd[i,j]= Sum[ (D[chrisddu[i,j,n],x[[n]]] - D[chrisddu[j,n,n],x[[i]]]), {n,1,2}] ]];


For[i=1,i<3,i++,
	For[j=i,j<3,j++,
		 ric2dd[j,i]=ric2dd[i,j]=
								Sum[ Sum[(chrisddu[k,n,n]*chrisddu[i,j,k] - chrisddu[i,k,n]*chrisddu[n,j,k]) , {k,1,2}], {n,1,2}] ]];

For[i=1,i<3,i++,
	For[j=i,j<3,j++,
			 ricdd[j,i]=ricdd[i,j]=ric1dd[i,j]+ric2dd[i,j]]];

gdd //MatrixForm
guu //MatrixForm

(* Define the vector v and scalar v0 required for the Einstein equations. See Chapter 3, equations 3.3 and 3.4. *)
vvecd=Array[vv,{2,1}];For[i=1,i<3,i++,vvecd[[i]]=D[Exp[p]*S[z],x[[i]]]/(Exp[p]*S[z])];
vvecd=Simplify[vvecd];

vvecu=Simplify[guu.vvecd];
V0=Simplify[-1/(Exp[2p]*S[z]^2)+2Sum[D[vvecu[[i]],x[[i]]]+Sum[chrisddu[i,j,i]vvecu[[j]],{j,1,2}],{i,1,2}]+3Sum[vvecd[[i]]vvecu[[i]],{i,1,2}]];
Clear[al, be, Ga, Gga, g, gg, ggg, ksc, kksc, kscc, all, bee, u, v, w, delrho, rho]

Clear[ ddelrho,delrho, rrho,rho, ksc, g, al, be, Ga, ga]

(* The background density (see Chapter 2, equation 2.4): *)
rho[z_, p_]=q[z]*Exp[-2*p];

(* The density perturbation (see Chapter 5, equation 5.1): *)
delrho[z_, p_]=Exp[-4*p]*ddelrho[z,p];

(* The scalar metric perturbation for even parity (see Chapter 3, equation 3.15):  *)
ksc[z_,p_]:=kksc[z,p]*Exp[-2*p];

(* The perturbation of the dust velocity (see Chapter 5, equation 5.2): *)
g[p_]=Exp[-p]*gg[p];
Ga[z_,p_]=Exp[-p]*Gga[z,p];

(* We define ga[z,p] in terms of Ga[z,p] and an initial data function g[p], see Chapter 5, equation 5.4. *)
ga[z_,p_]=Ga[z,p]+g[p];

(* The tensor stress-energy perturbation: *)
T00=2*rho[z,p]*Exp[p]*D[Ga[z,p], z]+Exp[2*p]*delrho[z,p];
T01=Exp[p]*rho[z,p]*(D[Ga[z,p], p]+z*D[Ga[z,p], z])+Exp[2*p]*z*delrho[z,p];
T11=2*rho[z,p]*z*Exp[p]*D[Ga[z,p], p]+delrho[z,p]*z^2*Exp[2*p];
delttensor={{T00, T01}, {T01, T11}};

(* The vector stress-energy perturbation: *)
v1=Exp[p]*rho[z,p]*ga[z,p];
v2=Exp[p]*rho[z,p]*ga[z,p]*z;
deltvector={v1, v2};

(* The scalar stress-energy perturbation: *)
T2=0;
T3=0;

Clear[de]

(* The tensor metric perturbation, see equation 3.15: *)
ktdd={{al[z,p],be[z,p]},{be[z,p],de[z,p]}};
ktdu=ktdd.guu;
ktuu=guu.ktdu;

(*Various derivatives of the tensor metric perturbation: *)
For[i=1,i<3,i++,divktd[i]=Sum[D[ktdu[[i,j]],x[[j]]]-Sum[chrisddu[i,j,k]*ktdu[[k,j]]-chrisddu[k,j,j]*ktdu[[i,k]],{k,1,2}],{j,1,2}]];
For[i=1,i<3,i++,divktu[i]=Sum[guu[[i,j]]divktd[j],{j,1,2}]];
trk=Sum[guu[[i,j]]ktdd[[i,j]],{i,1,2},{j,1,2}];
For[i=1,i<3,i++,gradtrku[i]=Sum[guu[[i,j]]D[trk,x[[j]]],{j,1,2}]];
boxtrk=Sum[D[gradtrku[i],x[[i]]]+Sum[chrisddu[i,j,i]gradtrku[j],{j,1,2}],{i,1,2}];
For[i=1,i<3,i++,gradkscu[i]=Sum[guu[[i,j]]D[ksc[z,p],x[[j]]],{j,1,2}]];
boxksc=Sum[D[gradkscu[i],x[[i]]]+Sum[chrisddu[i,j,i]gradkscu[j],{j,1,2}],{i,1,2}];

(*These terms should be the same: *)
test=Simplify[trk-Sum[ktdu[[i,i]],{i,1,2}]]

Clear[de]
eqn1=trk-2T2;
Simplify[eqn1];
Simplify[Solve[eqn1==0,de[z,p]]]

(* Eliminate de(z,p): *)
de[z_,p_]=2 z be[z,p]+al[z,p] (((S[z])^2)-2 z S[z] S'[z]+(z^2) (-1+((S'[z])^2)));

For[i=1,i<3,i++,eqn2[i]=Sum[D[ktdu[[i,j]],x[[j]]]-Sum[chrisddu[i,j,k]*ktdu[[k,j]]-chrisddu[k,j,j]*ktdu[[i,k]],{k,1,2}],{j,1,2}]-D[trk,x[[i]]]+trk*vvecd[[i]]-D[ksc[z,p],x[[i]]]-16*Pi*deltvector[[i]]];

(* The first component: *)
E1:=Simplify[eqn2[1]]
E1=Simplify[E1]*(((S[z]-z S'[z]))^3)*Exp[2*p]

(* The second component: *)
E2:=Simplify[eqn2[2]]
E2=E2*(((S[z]-z S'[z]))^3)*Exp[2*p]

eqn3=-Sum[D[divktu[i],x[[i]]]+Sum[chrisddu[i,j,i]divktu[j],{j,1,2}],{i,1,2}]+boxtrk-2Sum[divktu[i]vvecd[[i]],{i,1,2}]+Sum[gradtrku[i]vvecd[[i]],{i,1,2}]+Sum[ric4uu[i,j](ktdd[[i,j]]-ksc[z,p]gdd[[i,j]]),{i,1,2},{j,1,2}]-(L(L+1)/(2*S[z]^2*Exp[2*p]))trk+boxksc+2Sum[gradkscu[i]vvecd[[i]],{i,1,2}]-2T3;

E3:=Simplify[eqn3]
E3=E3*S[z]*(((S[z]-z S'[z]))^3)*Exp[4*p]

eqn3alt=-boxksc -2Sum[guu[[i,j]](D[deltvector[[i]],x[[j]]]-Sum[chrisddu[i,j,k]deltvector[[k]],{k,1,2}]),{i,1,2},{j,1,2}]+boxtrk-2Sum[divktu[i]vvecd[[i]],{i,1,2}]+Sum[gradtrku[i]vvecd[[i]],{i,1,2}]+Sum[ric4uu[i,j](ktdd[[i,j]]-ksc[z,p]gdd[[i,j]]),{i,1,2},{j,1,2}]-(L(L+1)/(2*S[z]^2*Exp[2*p]))trk+boxksc+2Sum[gradkscu[i]vvecd[[i]],{i,1,2}]-2T3;

E3alt:=Simplify[eqn3alt]
E3alt=E3alt*S[z]*(((S[z]-z S'[z]))^3)*(1/2)*Exp[4*p]

For[i=1,i<3,i++,For[j=i,j<3,j++,eqn4[i,j]=Sum[vvecu[[k]](D[ktdd[[k,i]],x[[j]]]+D[ktdd[[k,j]],x[[i]]]-D[ktdd[[i,j]],x[[k]]]-2Sum[chrisddu[i,j,l]ktdd[[k,l]],{l,1,2}]),{k,1,2}]-gdd[[i,j]]*Sum[vvecu[[k]](2divktd[k]-D[trk,x[[k]]]),{k,1,2}]-(D[ksc[z,p],x[[i]]]*vvecd[[j]]+D[ksc[z,p],x[[j]]]*vvecd[[i]]+D[ksc[z,p],x[[i]],x[[j]]]-Sum[chrisddu[i,j,k]D[ksc[z,p],x[[k]]],{k,1,2}])+ktdd[[i,j]]*(V0+(L(L+1))/(2*S[z]^2*Exp[2*p]))-gdd[[i,j]]*(L(L+1)trk/(2*S[z]^2*Exp[2*p])+2Sum[ktdu[[k,l]]*(D[vvecu[[k]],x[[l]]]+Sum[chrisddu[l,m,k]vvecu[[m]],{m,1,2}]),{k,1,2},{l,1,2}]+3Sum[ktdd[[k,l]]vvecu[[k]]vvecu[[l]],{k,1,2},{l,1,2}]-boxksc-3Sum[vvecu[[k]]D[ksc[z,p],x[[k]]],{k,1,2}]+(L-1)(L+2)ksc[z,p]/(2*S[z]^2*Exp[2*p]))-8*Pi*delttensor[[i,j]]]];

(* The first component: *)
E4:=Simplify[eqn4[1,1]]
E4=E4*2*S[z]*S[z] (((S[z]-z S'[z]))^3)*Exp[2*p]

Simplify[Coefficient[E4,ddelrho[z,p]]]

Simplify[Solve[E4==0,ddelrho[z,p]]]

ddelrho[z_,p_]=(1/(16 π ((S[z])^2) (((S[z]-z S'[z]))^3)))(al[z,p] ((-2+L+(L^2)) ((S[z])^3)+3 (2+L+(L^2)) (z^2) S[z] ((S'[z])^2)-(-2+L+(L^2)) (z^3) ((S'[z])^3)-z ((S[z])^2) (3 (2+L+(L^2)) S'[z]+4 z S''[z]))+kksc[z,p] ((-2+L+(L^2)) ((S[z])^3)+3 (2+L+(L^2)) (z^2) S[z] ((S'[z])^2)-(-2+L+(L^2)) (z^3) ((S'[z])^3)-z ((S[z])^2) (3 (2+L+(L^2)) S'[z]+4 z S''[z]))+2 (2 be[z,p] (-(z^2) ((S'[z])^3)+((S[z])^2) (S'[z]+z S''[z]))+S[z] (-16 π q[z] ((S[z])^4) Derivative[1,0][Gga][z,p]+((S[z])^3) (-z S''[z] Derivative[1,0][kksc][z,p]+S'[z] (Derivative[1,0][al][z,p]+48 π z q[z] Derivative[1,0][Gga][z,p]+Derivative[1,0][kksc][z,p]))-z ((S'[z])^2) (z Derivative[0,1][al][z,p]-2 Derivative[0,1][be][z,p]+z (3 Derivative[0,1][kksc][z,p]+z (-1+((S'[z])^2)) Derivative[1,0][al][z,p]+2 Derivative[1,0][be][z,p]-3 z Derivative[1,0][kksc][z,p]+z ((S'[z])^2) Derivative[1,0][kksc][z,p]))+((S[z])^2) (Derivative[0,1][al][z,p]+2 Derivative[0,1][kksc][z,p]-Derivative[0,2][kksc][z,p]-z Derivative[1,0][al][z,p]-3 z ((S'[z])^2) Derivative[1,0][al][z,p]-48 π (z^2) q[z] ((S'[z])^2) Derivative[1,0][Gga][z,p]-3 z Derivative[1,0][kksc][z,p]-3 z ((S'[z])^2) Derivative[1,0][kksc][z,p]+2 (z^2) S'[z] S''[z] Derivative[1,0][kksc][z,p]+2 z Derivative[1,1][kksc][z,p]-(z^2) Derivative[2,0][kksc][z,p])+S[z] (-(z^3) ((S'[z])^2) S''[z] Derivative[1,0][kksc][z,p]+(z^2) ((S'[z])^3) (3 Derivative[1,0][al][z,p]+16 π z q[z] Derivative[1,0][Gga][z,p]+3 Derivative[1,0][kksc][z,p])+(z^2) S''[z] (Derivative[0,1][kksc][z,p]-z Derivative[1,0][kksc][z,p])+S'[z] (-2 Derivative[0,1][be][z,p]+z (Derivative[0,1][kksc][z,p]+Derivative[0,2][kksc][z,p]+2 Derivative[1,0][be][z,p]-2 z Derivative[1,1][kksc][z,p]+(z^2) Derivative[2,0][kksc][z,p]))))));

(* The second component: *)
E5:=Simplify[eqn4[1,2]]
E5=E5*2 ((S[z])^2) (((S[z]-z S'[z]))^3)*Exp[2*p]

(* The third component: *)
E6:=Simplify[eqn4[2,2]]
E6=Simplify[E6*2*S[z]^2*(((S[z]-z S'[z]))^1)]*Exp[2*p]

Clear[EQ1, EQ2, EQ4, EQ5, EQ6, EQ7]

(* The initial system: *)
EQ1={E1,E2,E5,E6};

(* The first order reduction: *)EQ2=EQ1/.{kksc[z,p]->u[z,p],Derivative[1,0][kksc][z,p]->w[z,p],Derivative[0,1][kksc][z,p]->v[z,p],Derivative[1,1][kksc][z,p]->D[w[z,p],p],Derivative[2,0][kksc][z,p]->D[w[z,p],z],Derivative[0,2][kksc][z,p]->D[v[z,p],p]};

(* Add in the stress-energy conservation equation: *)
EQ4=Append[EQ2, D[Gga[z,p], z]+(1/2)*al[z,p]];

(* Add in the trivial evolution equations arising from the first order reduction: *)
EQ5=Append[EQ4,D[u[z,p],z]-w[z,p]];
EQ6=Append[EQ5,D[v[z,p],z]-D[w[z,p],p]];

(* Define the initial state vector X: *)
X={al[z,p],be[z,p],u[z,p],v[z,p],w[z,p], Gga[z,p]}

Dimensions[EQ6]
Dimensions[X]

A=Array[amat,{7,6}];For[i=1,i<8,i++,For[j=1,j<7,j++,A[[i,j]]=Coefficient[EQ6[[i]],D[X[[j]],z]]]];
B=Array[bmat,{7,6}];For[i=1,i<8,i++,For[j=1,j<7,j++,B[[i,j]]=Coefficient[EQ6[[i]],D[X[[j]],p]]]];
F=Array[cmat,{7,6}];For[i=1,i<8,i++,For[j=1,j<7,j++,F[[i,j]]=Coefficient[EQ6[[i]],X[[j]]]]];

Simplify[A] //MatrixForm

checkA2=(({{-((S[z])^3)+3 z ((S[z])^2) S'[z]+(z^2) S[z] (1-3 ((S'[z])^2))+(z^3) S'[z] (-1+((S'[z])^2)), z (-S[z]+z S'[z]), 0, 0, 0, 0}, {z (-((S[z])^3)+3 z ((S[z])^2) S'[z]+(z^2) S[z] (1-3 ((S'[z])^2))+(z^3) S'[z] (-1+((S'[z])^2))), -((S[z])^3)+3 z ((S[z])^2) S'[z]+(z^3) S'[z] (1+((S'[z])^2))-(z^2) S[z] (1+3 ((S'[z])^2)), 0, 0, 0, 0}, {2 ((S[z])^2) (((S[z]-z S'[z]))^3), 0, 0, 0, 2 z ((S[z])^2) (((S[z]-z S'[z]))^3), 16 π z q[z] ((S[z])^2) (((S[z]-z S'[z]))^3)}, {2 S[z] (-((S[z])^3) S'[z]-3 (z^2) S[z] ((S'[z])^3)+(z^3) ((S'[z])^2) (-1+((S'[z])^2))+((S[z])^2) (z+3 z ((S'[z])^2))), 4 S[z] (((S[z]-z S'[z]))^2), 0, 0, -2 ((S[z])^2) (((S[z])^3)-3 z ((S[z])^2) S'[z]-(z^3) S'[z] (-2+((S'[z])^2))+(z^2) S[z] (-2+3 ((S'[z])^2))), 32 π (z^2) q[z] ((S[z])^2) (S[z]-z S'[z])}, {0, 0, 0, 0, 0, 1}, {0, 0, 1, 0, 0, 0}}));

(* This matrix has vanishing determinant: *)
Simplify[Det[checkA2]]

(* The coefficients: *)
Clear[c1, c2, c3, c5, c6]
sum=c1*Part[checkA2, 1]+c2*Part[checkA2, 2]+c3*Part[checkA2, 3]-Part[checkA2, 4]+c5*Part[checkA2, 5]+c6*Part[checkA2, 6];
Solve[sum==0,{c1,c2,c3, c5, c6}]

Simplify[sum/.{c5->-16 (-π q[z] ((S[z])^5)+3 π z q[z] ((S[z])^4) S'[z]-3 π (z^2) q[z] ((S[z])^3) ((S'[z])^2)+π (z^3) q[z] ((S[z])^2) ((S'[z])^3)),c6->0,c1->((2 (-(z^2) S[z]+((S[z])^3)-2 z ((S[z])^2) S'[z]+(z^2) S[z] ((S'[z])^2)))/(z (-S[z]+z S'[z]))),c2->-((2 S[z])/(S[z]-z S'[z])),c3->-((-2 (z^2)+((S[z])^2)-2 z S[z] S'[z]+(z^2) ((S'[z])^2))/(z (((-S[z]+z S'[z]))^2)))}]

Clear[constraint]
constraint=Simplify[c1*EQ6[[1]]+c2*EQ6[[2]]+c3*EQ6[[3]]-EQ6[[4]]+c5*EQ6[[5]]+c6*EQ6[[6]]/.{c5->-16 (-π q[z] ((S[z])^5)+3 π z q[z] ((S[z])^4) S'[z]-3 π (z^2) q[z] ((S[z])^3) ((S'[z])^2)+π (z^3) q[z] ((S[z])^2) ((S'[z])^3)),c6->0,c1->((2 (-(z^2) S[z]+((S[z])^3)-2 z ((S[z])^2) S'[z]+(z^2) S[z] ((S'[z])^2)))/(z (-S[z]+z S'[z]))),c2->-((2 S[z])/(S[z]-z S'[z])),c3->-((-2 (z^2)+((S[z])^2)-2 z S[z] S'[z]+(z^2) ((S'[z])^2))/(z (((-S[z]+z S'[z]))^2)))}];
constraint=Simplify[constraint*z*(((S[z]-z S'[z]))^(-2)) ]

S[z_]=sform[z];
q[z_]=qform[z];
constraintz=Simplify[constraint]
Clear[S, q]

(* Here the Q_{i}'s are arbitrary coefficients depending on z only: *)
Clear[Q1, Q2, Q3, Q4, Q5, Q6, Q7, EQ7]
EQ7={EQ6[[1]]+Q1[z]*constraint, EQ6[[2]]+Q2[z]*constraint,EQ6[[3]]+Q3[z]*constraint, EQ6[[5]]+Q5[z]*constraint,
EQ6[[6]]+Q6[z]*constraint, EQ6[[7]]+Q7[z]*constraint};

(* We identify the matrix coefficients J, G and H in the system J(z) ∂X / ∂z + G(z) ∂X / ∂p + H(z) X = Source. This corresponds to equation 5.6 in Chapter 5.  *)
J=Array[jmat,{6,6}];For[i=1,i<7,i++,For[j=1,j<7,j++,J[[i,j]]=Coefficient[EQ7[[i]],D[X[[j]],z]]]];
G=Array[gmat,{6,6}];For[i=1,i<7,i++,For[j=1,j<7,j++,G[[i,j]]=Coefficient[EQ7[[i]],D[X[[j]],p]]]];
H=Array[hmat,{6,6}];For[i=1,i<7,i++,For[j=1,j<7,j++,H[[i,j]]=Coefficient[EQ7[[i]],X[[j]]]]];

Simplify[J] //MatrixForm

Simplify[Det[J]]

(* Identify the source vector: *)
Source=Simplify[EQ7-J.D[X,z]-G.D[X,p]-H.X]

(* Write the system in a slightly different form: *)
Xdot=-Simplify[Inverse[J].(G.D[X,p]+H.X+Source)];
Dimensions[Xdot]

(*The constraint should be propagated for any choice of functions Q_{i}(z). However, to simplify matters, we make the simplest possible choice of Q_{i}'s: *)
Q1[z_]=0;
Q2[z_]=0;
Q3[z_]=0;
Q4[z_]=0;
Q5[z_]=0;
Q6[z_]=0;
Q7[z_]=0;

condot=Simplify[D[constraint,z]];

(* Insert the forms of the X derivatives into the derivative of the constraint: *)
condot1=condot/.{D[X[[1]],z]->Xdot[[1]],D[X[[2]],z]->Xdot[[2]],D[X[[3]],z]->Xdot[[3]],D[X[[4]],z]->Xdot[[4]],D[X[[5]],z]->Xdot[[5]], D[X[[6]],z]->Xdot[[6]]};

(* Insert all second order X derivatives: *)
condot2=condot1/.{D[X[[1]],z,p]->D[Xdot[[1]],p],D[X[[2]],z,p]->D[Xdot[[2]],p],D[X[[3]],z,p]->D[Xdot[[3]],p],D[X[[4]],z,p]->D[Xdot[[4]],p],D[X[[5]],z,p]->D[Xdot[[5]],p], D[X[[6]],z,p]->D[Xdot[[6]],p]};

conprime=Simplify[D[constraint,p]];

conprime

Clear[t1, testevolvecon, testevolvecon2]
t1=Simplify[Coefficient[conprime,Derivative[0,2][Gga][z,p]]];
t2=Simplify[Coefficient[condot2, Derivative[0,2][Gga][z,p]]];

testevolvecon=t1*condot2-t2*conprime;

t3=Simplify[Coefficient[constraint, Derivative[0,1][Gga][z,p]]];
t4=Simplify[Coefficient[testevolvecon, Derivative[0,1][Gga][z,p]]];

S[z_]=sform[z];
q[z_]=qform[z];

testevolvecon2=Simplify[t3*testevolvecon-t4*constraint]

(* Insert the trivial constraint: *)
testevolvecon3=testevolvecon2/.{Derivative[0,1][u][z,p]->v[z,p]}

Clear[S, q]

(* By multiplying through by J^{-1}, we can write the system as ∂ X/∂ z + Q ∂ X/∂ p + T X =  SS*)
Q=Simplify[Inverse[J].G];
T=Simplify[Inverse[J].H];
SS=Simplify[Inverse[J].Source];

SS=Simplify[Inverse[J].Source];

SS //MatrixForm

Simplify[Q] //MatrixForm

(* sim is the similarity matrix for Q and jordQ is the Jordan canonical form of Q. *)
{sim, jordQ}=Simplify[JordanDecomposition[Q]];

Simplify[jordQ] //MatrixForm

Eigenvalues[jordQ]

siminv=Inverse[sim];

V=Simplify[siminv.(T.sim+D[sim,z])];

xx={al[z,p],be[z,p],kscc[z,p],Derivative[0,1][kscc][z,p],Derivative[1,0][kscc][z,p],Gga[z,p]}

(* Define Y: *)
Y=Simplify[Inverse[sim].X]

Clear[y1, y2, y3, y4, y5, y6, all, bee, u, v, w, Gga]
y={y1[z,p],y2[z,p],y3[z,p],y4[z,p],y5[z,p],y6[z,p]};

(* We solve for the components of X in terms of the components of Y: *)
Simplify[Solve[{y==Y},X]]

(* Replace the components of X with the components of Y: *)
v[z_,p_]=y1[z,p]+z y4[z,p]+z y5[z,p]-S[z] y5[z,p]+z y6[z,p]+S[z] y6[z,p]+z y5[z,p] S'[z]-z y6[z,p] S'[z];
u[z_,p_]=y3[z,p];
w[z_,p_]=-y2[z,p]+y4[z,p]+y5[z,p]+y6[z,p];
Gga[z_,p_]=((y2[z,p])/(8 π q[z]));
al[z_,p_]=-((S[z] (y5[z,p] (-1+S'[z])+y6[z,p] (1+S'[z])))/(-1+((S'[z])^2)));be[z_,p_]=-((S[z] (y5[z,p] (-1+S'[z]) (z-S[z]+z S'[z])+y6[z,p] (z-z ((S'[z])^2)+S[z] (1+S'[z]))))/(-1+((S'[z])^2)));

(* Rewrite the Einstein constraint in terms of Y. This corresponds to equation 5.8 in Chapter 5. *)
constrainty=Simplify[constraint]

(* Rewrite the trivial constraint in terms of Y. This corresponds to equation 5.9 in Chapter 5. *)
trivconstraint=Simplify[v[z,p]-Derivative[0,1][u][z,p]]

(* The equation of motion for Y is: *)
newSource=Simplify[siminv.SS];
ydot = - jordQ.D[y,p]-V.y-newSource;

Simplify[Solve[trivconstraint==0, y1[z,p]]]

y1[z_,p_]=-z y4[z,p]-z y6[z,p]-S[z] y6[z,p]+z y6[z,p] S'[z]-y5[z,p] (z-S[z]+z S'[z])+Derivative[0,1][y3][z,p];

(* We now reduce the system to a five dimensional system, with one remaining (non-trivial or Einstein) constraint. *)
Simplify[ydot]

(* The 5d system: *)
Clear[ww1, ww2, ww3, ww4, ww5]
wwdot={ydot[[2]], ydot[[3]], ydot[[4]], ydot[[5]], ydot[[6]]};
WW={ww1[z,p],ww2[z,p], ww3[z,p], ww4[z,p], ww5[z,p]};
y2[z_,p_]=ww1[z,p];
y3[z_,p_]=ww2[z,p];
y4[z_,p_]=ww3[z,p];
y5[z_,p_]=ww4[z,p];
y6[z_,p_]=ww5[z,p];
Simplify[wwdot]

PPP=Array[pmat,{5,5}];For[i=1,i<6,i++,For[j=1,j<6,j++,PPP[[i,j]]=Coefficient[wwdot[[i]],D[WW[[j]],p]]]];
PPP //MatrixForm

Simplify[PPP] //MatrixForm

{PPPsim, PPPjord}=Simplify[JordanDecomposition[PPP]]

(* Diagonalize the system using the above similarity matrix.  *)
M=Simplify[Inverse[PPPsim].WW]

(* Solve for the components of W in terms of the components of M. *)
Clear[mm1, mm2, mm3, mm4, mm5]
mm={mm1[z,p], mm2[z,p], mm3[z,p], mm4[z,p], mm5[z,p]}

Simplify[Solve[{mm==M}, WW]]

ww3[z_,p_]=((mm3[z,p] (-1+S'[z]) (S[z]-z S'[z])+2 mm1[z,p] (-z ((S'[z])^2)+S[z] (S'[z]+z S''[z])))/((-1+S'[z]) (S[z]-z S'[z])));
ww1[z_,p_]=mm2[z,p];
ww4[z_,p_]=((mm4[z,p] (-1+S'[z])+mm1[z,p] (1+S'[z]))/(-1+S'[z]));
ww5[z_,p_]=mm1[z,p]+mm5[z,p];ww2[z_,p_]=-((2 mm1[z,p] S[z])/(-1+S'[z]));

constraint5d=Simplify[constrainty]

Clear[S,q]
Clear[mm1, mm2, mm3, mm4, mm5]
mm

(* We now identify the matrix coefficients and source term in MJ(z) ∂ M / ∂ z + MG(z) ∂ M / ∂ p + MH(z) M = msource *)
mmeq={D[ww1[z,p],z]-wwdot[[1]],D[ww2[z,p],z]-wwdot[[2]],D[ww3[z,p],z]-wwdot[[3]], D[ww4[z,p],z]-wwdot[[4]],D[ww5[z,p],z]-wwdot[[5]] };

MJ=Array[mjmat,{5,5}];For[i=1,i<6,i++,For[j=1,j<6,j++,MJ[[i,j]]=Coefficient[mmeq[[i]],D[mm[[j]],z]]]];
MG=Array[mgmat,{5,5}];For[i=1,i<6,i++,For[j=1,j<6,j++,MG[[i,j]]=Coefficient[mmeq[[i]],D[mm[[j]],p]]]];
MH=Array[mhmat,{5,5}];For[i=1,i<6,i++,For[j=1,j<6,j++,MH[[i,j]]=Coefficient[mmeq[[i]],mm[[j]]]]];

msource={Simplify[Coefficient[mmeq[[1]], gg[p]]]*gg[p], Simplify[Coefficient[mmeq[[2]], gg[p]]]*gg[p], Simplify[Coefficient[mmeq[[3]], gg[p]]]*gg[p], Simplify[Coefficient[mmeq[[4]], gg[p]]]*gg[p], Simplify[Coefficient[mmeq[[5]], gg[p]]]*gg[p]};

msource

Simplify[MatrixForm[MJ]]

Simplify[MatrixForm[MG]]

Simplify[Det[MJ]]

S[z_]=sform[z];
Simplify[Det[MJ]]
Clear[S]

MGG=Inverse[MJ].MG;
MHH=Inverse[MJ].MH;

(* Note the fact that this matrix is symmetric which implies that this version of the system is symmetric hyperbolic. This is the matrix \tilde{A} which appears in equation 5.10, the form of the system which is used in Chapters 5 and 6.  \tilde{A} is listed below equation 5.10 in Chapter 5.  *)
Simplify[MGG] //MatrixForm

(* This is the matrix \tilde{C} which appears in equation 5.10 and in Chapters 5 and 6. Its components are given in Appendix B1. *)
Simplify[MHH] //MatrixForm

(* This is the source term tilde{\Sigma} which appears in equation 5.10. Its components are given in Appendix B1. *)
newmsource=Simplify[Inverse[MJ].msource]

(* Check that everything is working - the Jordan decomposition of MGG should be MGG again. *)
{msimmer,mjord}=JordanDecomposition[MGG];
Simplify[mjord - MGG]

Clear[mm2, mm1, mm3, mm4, mm5, S]
Simplify[Solve[constraint5d==0, {mm2[z,p]}]]

mm2[z_,p_]=(mm4[z,p] (-1+S'[z]) (-S[z]+z S'[z]) (-(z^2) S'[z] (1+S'[z]) (-2+L+(L^2)-2 S'[z]-2 ((S'[z])^2))+((S[z])^3) (-8 π z q[z]+4 S''[z])+z S[z] (-2+L+(L^2)-4 ((S'[z])^3)-2 z S''[z]+2 S'[z] (-3+L+(L^2)+z S''[z])+((S'[z])^2) (-6+4 z S''[z]))-((S[z])^2) (-2+L+(L^2)-2 ((S'[z])^2)+4 z S''[z]+S'[z] (-2-8 π (z^2) q[z]+8 z S''[z])))+(1+S'[z]) (mm5[z,p] (-S[z]+z S'[z]) ((z^2) (-1+S'[z]) S'[z] (-2+L+(L^2)+2 S'[z]-2 ((S'[z])^2))-4 ((S[z])^3) (2 π z q[z]+S''[z])+z S[z] (-2+L+(L^2)+4 ((S'[z])^3)+2 z S''[z]-2 S'[z] (-3+L+(L^2)-z S''[z])-2 ((S'[z])^2) (3+2 z S''[z]))+((S[z])^2) (-2+L+(L^2)-2 ((S'[z])^2)-4 z S''[z]+S'[z] (2+8 π (z^2) q[z]+8 z S''[z])))+2 (8 π z mm1[z,p] q[z] ((S[z])^4)-4 mm1[z,p] ((S[z])^3) S'[z]-16 π (z^2) mm1[z,p] q[z] ((S[z])^3) S'[z]+12 z mm1[z,p] ((S[z])^2) ((S'[z])^2)+8 π (z^3) mm1[z,p] q[z] ((S[z])^2) ((S'[z])^2)-12 (z^2) mm1[z,p] S[z] ((S'[z])^3)+4 (z^3) mm1[z,p] ((S'[z])^4)+16 π gg[p] q[z] (-1+S'[z]) (((-S[z]+z S'[z]))^3)+2 z mm1[z,p] ((S[z])^3) S''[z]-2 (z^2) mm1[z,p] ((S[z])^2) S'[z] S''[z]+2 (z^3) mm1[z,p] ((S[z])^2) ((S''[z])^2)+mm3[z,p] S[z] (-1+S'[z]) (S[z]-z S'[z]) (S[z]+z (-S'[z]+z S''[z]))+((S[z])^3) Derivative[0,1][mm3][z,p]-2 z ((S[z])^2) S'[z] Derivative[0,1][mm3][z,p]-((S[z])^3) S'[z] Derivative[0,1][mm3][z,p]+(z^2) S[z] ((S'[z])^2) Derivative[0,1][mm3][z,p]+2 z ((S[z])^2) ((S'[z])^2) Derivative[0,1][mm3][z,p]-(z^2) S[z] ((S'[z])^3) Derivative[0,1][mm3][z,p])))/(2 (S[z]-z S'[z]) (-1+((S'[z])^2)) (2 ((S[z])^2)+2 (z^2) ((S'[z])^2)+z S[z] (-4 S'[z]+z S''[z])));

(*The m equation: *)
mdot=-(MGG.D[mm,p]+MHH.mm);

(* Construct the source term here  *)
newmsource

(* The 4-dimensional source: *)
mmsource4d={0,((16 π gg[p] q[z] (-S[z]+z S'[z]))/(z S[z])),((8 π gg[p] q[z] (1+S'[z]) (S[z]-z S'[z]))/(S[z] (z-S[z]+z S'[z]))),-((8 π gg[p] q[z] (-1+S'[z]) (S[z]-z S'[z]))/(S[z] (z+S[z]-z S'[z])))};

(* The four-dimensional system: *)
nneq={D[mm1[z,p],z]-mdot[[1]],D[mm3[z,p],z]- mdot[[3]], D[mm4[z,p],z]-mdot[[4]],D[mm5[z,p],z]- mdot[[5]]};
NN={n1[z,p], n2[z,p], n3[z,p], n4[z,p]};
mm1[z_, p_]=n1[z,p];
mm3[z_, p_]=n2[z,p];
mm4[z_, p_]=n3[z,p];
mm5[z_, p_]=n4[z,p];

Clear[S, q]
Qmatnoz=Array[qmat,{4,4}];For[i=1,i<5,i++,For[j=1,j<5,j++,Qmatnoz[[i,j]]=Coefficient[nneq[[i]],D[NN[[j]],z]]]];
Pmatnoz=Array[ppmat,{4,4}];For[i=1,i<5,i++,For[j=1,j<5,j++,Pmatnoz[[i,j]]=Coefficient[nneq[[i]],D[NN[[j]],p]]]];
Tmatnoz=Array[ttmat,{4,4}];For[i=1,i<5,i++,For[j=1,j<5,j++,Tmatnoz[[i,j]]=Coefficient[nneq[[i]],NN[[j]]]]];

S[z_]=sform[z];
q[z_]=qform[z];

(* The Q, P and T matrices with background information included: *)
Qmat=Array[qmat,{4,4}];For[i=1,i<5,i++,For[j=1,j<5,j++,Qmat[[i,j]]=Coefficient[nneq[[i]],D[NN[[j]],z]]]];
Pmat=Array[pmat,{4,4}];For[i=1,i<5,i++,For[j=1,j<5,j++,Pmat[[i,j]]=Coefficient[nneq[[i]],D[NN[[j]],p]]]];
Pmat=Simplify[Pmat] ;
Tmat=Array[tmat,{4,4}];For[i=1,i<5,i++,For[j=1,j<5,j++,Tmat[[i,j]]=Coefficient[nneq[[i]],NN[[j]]]]];

(* We can confirm that adding in z values makes no difference yet - but it will at the Jordan decomposition. *)
Simplify[Qmatnoz-Qmat]
Simplify[Pmatnoz-Pmat]
Simplify[Tmatnoz-Tmat]

(* We now multiply through by (Q)^{-1} to get the equation ∂ N/∂ z + newP ∂ N/∂ p + newT N = newnsource. *)
Clear[S, q]
newP=Inverse[Qmat].Pmat;
newPnoz=Inverse[Qmatnoz].Pmatnoz;

(* Calculate the Jordan decomposition again: *)
{newPsim,newPjord}=Simplify[JordanDecomposition[newP]]
{newPsimnoz, newPjordnoz}=Simplify[JordanDecomposition[newPnoz]]

(* Here is where the background information matters: newPjord and newPjordnoz are not the same. *)
S[z_]=sform[z];
q[z_]=qform[z];
Simplify[newPjord-newPnoz]

(* This matrix corresponds to the matrix E(z) appearing in equation 5.12 in Chapter 5. It is not symmetric, and so the system in the form 5.12 is not symmetric hyperbolic. The components of this matrix are listed in Appendix B1.*)
newPjord //MatrixForm

pp={pp1[z,p], pp2[z,p], pp3[z,p], pp4[z,p]};
PP=Simplify[Inverse[newPsim]. NN];
Simplify[Solve[{PP==pp}, NN]]