Cartesian Kerr-Schild Cartesian Tetrad

CKS Metric

ClearAll["Global`*"]
$Assumptions = _ ∈ Reals;

Mathematica does not begin indexing at 0; I find it most straightforward to produce fake indices for looping over metrics, four vectors,  spatial metrics, and three vectors that are consistent with the usual notation (e.g., the time component of a four vector being the 0 index).  In the following let f0-f3 be the loop limits for four row and/or column data structures and t1-t3 be the loop limits for three row and/or column data structures.

f0=1; f1=2; f2=3; f3=4;
t1=1; t2=2; t3=3;

We begin the with covariant form of the Cartesian Kerr-Schild metric (Visser 2008), 

(Subscript[g, μν])=(Subscript[η, μν])+f (Subscript[l, μ])(Subscript[l, ν])

where (Subscript[η, μν])is the Minkowski metric, and

f=((2 m (r^3))/((r^4)+(a^2)(z^2))),

and (Subscript[l, μ]) is a null vector with respect to (Subscript[g, μν])and (Subscript[η, μν])

(Subscript[l, μ])= (1, ((r x + ay)/((r^2)+(a^2))), ((r y-a x)/((r^2)+ (a^2))),(z/r))

and
r(x,y,z)=(Sqrt[(((R^2)-(a^2)+(Sqrt[((((R^2)-(a^2)))^2)+4(a^2)(z^2)]))/2)]),

where (R^2)=(x^2)+(y^2)+(z^2) and m and a are the mass and spin of the black hole respectively.  We also have

(l^μ)= (-1, ((r x + ay)/((r^2)+(a^2))), ((r y-a x)/((r^2)+ (a^2))),(z/r)).

By definition, 

 -1+((Subscript[l, 1])^2)+((Subscript[l, 2])^2)+((Subscript[l, 3])^2)=0.

ll={ll0,l1,l2,l3};
lu={lu0,l1,l2,l3};
η={{-1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
gcov=FullSimplify@Table[ f*ll[[μ]]*ll[[ν]],{μ,f0,f3},{ν,f0,f3}]+η /.{ll0-> 1};
gcon=FullSimplify@Table[-f*lu[[μ]]*lu[[ν]],{μ,f0,f3},{ν,f0,f3}]+η /.{lu0->-1};

Let’s confirm a few things about the metric.  First, let’s check that the inverse of (Subscript[g, μν])gives (g^(μν))

FullSimplify[Inverse[gcov]==gcon//.{(l1^2)->1-(l2^2)-(l3^2)}]

Second, let’s check that the (Sqrt[-g])=1

FullSimplify[Sqrt[-Det[gcov]]//.{(l1^2)->1-(l2^2)-(l3^2)}]

Equations Solved

Next, I’ll recap @cwhite’s notes on the equations solved in Athena++.  For now, we will restrict ourselves to vacuum gray transfer.  The equations solved (in discretized form) are

(Subscript[(1/(Δt)) [(Subscript[⟨(n^0)(Subscript[n, 0]) I ⟩, V Ω])] , (Subscript[x, -]^0)]^(Subscript[x, +]^0))   +  (∑_( i∈{1,2,3})) (Subscript[(1/(ΔV))[Δ (Subscript[A, i])(Subscript[⟨(n^i)(Subscript[n, 0]) I ⟩, t (Subscript[A, i]) Ω])], (Subscript[x, -]^i)]^(Subscript[x, +]^i)) + (∑_( a∈{ζ,ψ})) (Subscript[(1/(ΔΩ))[Δ (Subscript[λ, a])(Subscript[⟨(n^a)(Subscript[n, 0]) I ⟩, t V (Subscript[λ, a])])], (SubMinus[a])]^(SubPlus[a])) = 0

where Δt is the numerical timestep, I is the frequency integrated intensity, (n^μ) and (n^a) are normal vector coordinate coordinate components for the spatial and angular meshes, the second term corresponds to spatial fluxes on the Cartesian mesh, and the final term corresponds to angular fluxes on the angular mesh defined for each numerical cell.  Spatial fluxes are dependent on areas at cell faces (Subscript[ΔA, i]) and cell volumes ΔV.  Angular fluxes are dependent on arc lengths (Subscript[Δλ, a]) and solid angles ΔΩ.  

The infrastructure to compute arc lengths, solid angles, face areas, and cell volumes are already implemented in @cwhite’s version of GR radiation in Athena++.  However, we must make some extensions to compute (n^μ) and (n^a) when working with CKS. 

To move forward, we describe an implementation of an orthonormal tetrad for CKS.

Tetrad

(n^μ)  and (Subscript[n, 0]) are computed following

(n^μ) = (Subscript[e, ((OverHat[α]))]^μ) (n^(OverHat[α])) 
(Subscript[n, 0])  = (Subscript[g, 0μ ])(n^μ)

where (Subscript[e, ((OverHat[α]))]^μ) are orthonormal tetrad basis components and (n^(OverHat[α])) is defined in Davis & Gammie (2019) as

(n^(OverHat[α])) = (Subscript[δ, 0]^α) + sinζ cos ψ (Subscript[δ, 1]^α) + sinζ sin ψ (Subscript[δ, 2]^α) + cosζ  (Subscript[δ, 3]^α).

@cwhite’s implementation computes (n^a) as

(n^ζ)= (1/(sin ζ)) (n^(OverHat[α]))(n^(OverHat[β])) ((n^(OverHat[0])) (Subscript[ω, (OverHat[α]) (OverHat[β])]^(OverHat[ 3])) - (n^(OverHat[3])) (Subscript[ω, (OverHat[α]) (OverHat[β])]^(OverHat[ 0]))).   (**Note:** there appears to be some discrepancy between this and Equation 8 of DG2019, particularly by an additional entry of (n^(OverHat[0]))).  

(n^ζ)=(1/(((sin )^2)ζ)) (n^(OverHat[α]))(n^(OverHat[β])) ((n^(OverHat[2])) (Subscript[ω, (OverHat[α]) (OverHat[β])]^(OverHat[ 1])) - (n^(OverHat[1])) (Subscript[ω, (OverHat[α]) (OverHat[β])]^(OverHat[ 2]))). 

where 

(Subscript[ω, (OverHat[α]) (OverHat[β])]^(OverHat[ γ])) = (Subscript[e, ρ]^(((OverHat[γ])))) (Subscript[e, ((OverHat[β]))]^σ)  ((∂_(σ ))(Subscript[e, ((OverHat[α]))]^ρ) + (Subscript[Γ, σμ]^( ρ)) (Subscript[e, ((OverHat[α]))]^μ)),

(Subscript[e, ρ]^(((OverHat[α])))) = (η^(αβ)) (Subscript[g, ρσ ])(Subscript[e, ((OverHat[β]))]^σ) ,

(Subscript[Γ, ρσ]^( μ)) = (1/2)(g^(μλ))((∂_ρ) (Subscript[g, λσ]) + (∂_σ) (Subscript[g, λρ]) - (∂_λ) (Subscript[g, ρσ]))

and ∂ g and ∂ e correspond to derivatives of the metric and tetrad components, respectively.   Given an angular patch specified  by ζ and ψ, we can compute  (n^(OverHat[α])), and therefore , after specifying  (Subscript[e, ((OverHat[α]))]^μ) we can readily compute  (n^μ) and (n^a).   How should we assign (Subscript[e, ((OverHat[α]))]^μ)?

Again, following @cwhite’s notes, we seek orthonormal tetrad basis components (Subscript[e, ((OverHat[α]))]^μ) satisfying

(Subscript[g, μν ])(Subscript[e, ((OverHat[α]))]^μ) (Subscript[e, ((OverHat[β]))]^ν)= (Subscript[η, αβ]) 

We require that the timelike direction in the tetrad should be orthogonal to surfaces of constant time coordinate (x^0) and this direction should be future pointing:

(Subscript[e, ((OverHat[0]))]^μ)=-α (g^(0μ))

where α = (((-(g^00)))^(-1/2)).    The remaining components are computed via the Gram-Schmidt procedure.  We will consider a “Cartesian” tetrad obtained via Gram-Schmidt on vectors orthogonal to surfaces of constant (t,z,x,y), in that order.  The Gram-Schmidt procedure follows,

Let the spatial 3-metric be defined as 

(Subscript[γ, ij]) = (Subscript[g, ij])
(γ^ij)=(α^2)((g^(0i))(g^(0j))-(g^00) (g^ij)).  

Select coordinates (q^a) that are a function of the spatial coordinates (x^μ).   For a “Cartesian” tetrad, select (q^a)=(x,y,z) and (x^μ)=(t,x,y,z).  Define

(Subscript[A, ((OverHat[α]))i])= ((∂ (q^a))/(d (x^i)))

Calculate 

(Subscript[A, ((OverHat[3]))]^i)=(γ^ij)(Subscript[A, ((OverHat[3]))j])
(Subscript[B, ((OverHat[3]))]^i)=(Subscript[A, ((OverHat[3]))]^i)
(Subscript[C, ((OverHat[3]))]^i)= ((Subscript[γ, jk]) (Subscript[B, ((OverHat[3]))]^j)(Subscript[B, ((OverHat[3]))]^k))^(-1/2)(Subscript[B, ((OverHat[3]))]^i)

(Subscript[A, ((OverHat[1]))]^i)=(γ^ij)(Subscript[A, ((OverHat[1]))j])
(Subscript[B, ((OverHat[1]))]^i)=(Subscript[A, ((OverHat[1]))]^i)- (Subscript[γ, jk])(Subscript[A, ((OverHat[1]))]^j)(Subscript[C, ((OverHat[3]))]^k)(Subscript[C, ((OverHat[3]))]^i)
(Subscript[C, ((OverHat[1]))]^i)= ((Subscript[γ, jk]) (Subscript[B, ((OverHat[1]))]^j)(Subscript[B, ((OverHat[1]))]^k))^(-1/2)(Subscript[B, ((OverHat[1]))]^i)

(Subscript[A, ((OverHat[2]))]^i)=(γ^ij)(Subscript[A, ((OverHat[2]))j])
(Subscript[B, ((OverHat[2]))]^i)=(Subscript[A, ((OverHat[2]))]^i)- (Subscript[γ, jk])(Subscript[A, ((OverHat[2]))]^j)(Subscript[C, ((OverHat[3]))]^k)(Subscript[C, ((OverHat[3]))]^i)- (Subscript[γ, jk])(Subscript[A, ((OverHat[2]))]^j)(Subscript[C, ((OverHat[1]))]^k)(Subscript[C, ((OverHat[1]))]^i)
(Subscript[C, ((OverHat[1]))]^i)= ((Subscript[γ, jk]) (Subscript[B, ((OverHat[2]))]^j)(Subscript[B, ((OverHat[2]))]^k))^(-1/2)(Subscript[B, ((OverHat[2]))]^i)

Then set the orthonormal tetrad components following...

(Subscript[e, ((OverHat[a]))]^0)= 0
(Subscript[e, ((OverHat[a]))]^i)= (Subscript[C, ((OverHat[a]))]^i)

Gram-Schmidt

α=(((-gcon[[f0,f0]]))^(-1/2));
γcov=Table[gcov[[i,j]],{i,f1,f3},{j,f1,f3}];
γcon=(α^2)Table[(gcon[[f0,i]]*gcon[[f0,j]]-gcon[[f0,f0]]*gcon[[i,j]]),{i,f1,f3},{j,f1,f3}];

xμ={x,y,z};
q={x,y,z};

A3[idx_]:=Table[D[q[[t3]],xμ[[n]]],{n,t1,t3}][[idx]];  (*(Subscript[A, (3)idx]) = (dq^3)/(dx^idx)*)
matA3=FullSimplify@Table[Sum[γcon[[i,j]]*A3[j],{j,t1,t3}],{i,t1,t3}];
matB3=FullSimplify@matA3;
matC3=FullSimplify@Table[((Sum[γcov[[j,k]]*matB3[[j]]*matB3[[k]],{j,t1,t3},{k,t1,t3}])^(-1/2))matB3[[i]],{i,t1,t3}];
matC3=Assuming[f>0,FullSimplify[FullSimplify[matC3]/.{(-1+(l1^2)+(l2^2)+(l3^2))->0,((l1^2)+(l2^2)+(l3^2))->1}]];

A1[idx_]:=Table[D[q[[t1]],xμ[[n]]],{n,t1,t3}][[idx]];   (*(Subscript[A, (1)idx]) = (dq^1)/(dx^idx)*)
matA1=FullSimplify@Table[Sum[γcon[[i,j]]*A1[j],{j,t1,t3}],{i,t1,t3}];
matB1=FullSimplify@Table[matA1[[i]]-Sum[γcov[[j,k]]*matA1[[j]]*matC3[[k]]*matC3[[i]],{j,t1,t3},{k,t1,t3}],{i,t1,t3}];
matB1=Assuming[f>0,FullSimplify@matB1]//.{(-1+(l1^2)+(l2^2)+(l3^2))->0,(-1+(l1^2)+(l3^2))->-(l2^2)};
matB1=FullSimplify[matB1/.{(l2^2)->-(-1+(l1^2)+(l3^2))}];
matC1=Table[((Sum[γcov[[j,k]]*matB1[[j]]*matB1[[k]],{j,t1,t3},{k,t1,t3}])^(-1/2))matB1[[i]],{i,t1,t3}];
matC1=FullSimplify[matC1/.{(l3^2)->-(-1+(l1^2)+(l2^2))}];

A2[idx_]:=Table[D[q[[t2]],xμ[[n]]],{n,t1,t3}][[idx]];  (*(Subscript[A, (2)idx]) = (dq^2)/(dx^idx)*)
matA2=FullSimplify@Table[Sum[γcon[[i,j]]*A2[j],{j,t1,t3}],{i,t1,t3}];
matB2=Table[(matA2[[i]]
             -Sum[γcov[[j,k]]*matA2[[j]]*matC3[[k]]*matC3[[i]],{j,t1,t3},{k,t1,t3}]
             -Sum[γcov[[j,k]]*matA2[[j]]*matC1[[k]]*matC1[[i]],{j,t1,t3},{k,t1,t3}]),{i,t1,t3}];
matB2=Assuming[f>0,FullSimplify@matB2]/.(-1+(l1^2)+(l2^2)+(l3^2))->0;
matB2=FullSimplify[FullSimplify[matB2/.(l1^2)+(l2^2)->1-(l3^2)]/.{(-1+(l1^2)+(l2^2)+(l3^2))-> 0}];
matC2=Table[((Sum[γcov[[j,k]]*matB2[[j]]*matB2[[k]],{j,t1,t3},{k,t1,t3}])^(-1/2))matB2[[i]],{i,t1,t3}];
matC2=matC2/.{(l1^2)+(l2^2)->1-(l3^2)};

matC0=Table[-α gcon[[f0,μ]],{μ,f0,f3}];

Let’s clean up (Subscript[C, ((OverHat[a]))]^i) by eye, because Mathematica is not simplifying like I want it to...

matC0alt={(Sqrt[1+f]),-((f l1)/(Sqrt[1+f])),-((f l2)/(Sqrt[1+f])),-((f l3)/(Sqrt[1+f]))};
matC1alt={((Sqrt[1+f (l2^2)])/(Sqrt[1+f ((l1^2)+(l2^2))])),-((f l1 l2 )/((Sqrt[1+f (l2^2)])((Sqrt[1+f ((l1^2)+(l2^2))])))),0};
matC2alt={0,(1/(Sqrt[1+f (l2^2)])),0};
matC3alt={-((f l1 l3)/((Sqrt[1+f]) (Sqrt[1+f ((l1^2)+(l2^2))]))),-((f l2 l3)/((Sqrt[1+f]) (Sqrt[1+f ((l1^2)+(l2^2))]))),((Sqrt[1+f ((l1^2)+(l2^2))])/(Sqrt[1+f]))};

...and confirm that they are equivalent to the earlier iterations...

FullSimplify[matC0==matC0alt]
Assuming[f>0,FullSimplify[matC1==matC1alt]]
Assuming[f>0,FullSimplify[matC2==matC2alt]]
Assuming[f>0,FullSimplify[FullSimplify[matC3/.{(l3^2)->-(-1+(l1^2)+(l2^2))}]==matC3alt]]

tetrad={matC0alt,
        Prepend[matC1alt,0],
        Prepend[matC2alt,0],
        Prepend[matC3alt,0]};

Let’s cleanup even further by making the following substitutions...

{wa->(((1+f))^(1/2)),iwa->1/wa}
{wb->(1+f ((l1^2)+(l2^2)))^(1/2),iwb->1/wb}
{wc->(1+f (l2^2))^(1/2),iwc->1/wc}

tetrad=tetrad/.{(1+f)^(1/2)->wa,(1+f)^(-1/2)->iwa,
               (1+f ((l1^2)+(l2^2)))^(1/2)->wb,(1+f ((l1^2)+(l2^2)))^(-1/2)->iwb,
               (1+f (l2^2))^(1/2)->wc,(1+f (l2^2))^(-1/2)->iwc};

And at last, here is the Cartesian tetrad for CKS:

MatrixForm@tetrad

In the limit of M → 0, the tetrad should approach the identity matrix.

MatrixForm[
Limit[Evaluate[tetrad//.{wa->(((1+f))^(1/2)),iwa->1/wa,
                         wb->(1+f ((l1^2)+(l2^2)))^(1/2),iwb->1/wb,
                         wc->(1+f (l2^2))^(1/2),iwc->1/wc,
                         f->2 M (r^3)/((r^4)+(a^2)(z^2))}],{M->0}]]

Recall that we also need derivatives of the tetrad.  Let’s compute those here:

dwadx1->D[(((1+f[x]))^(1/2)),x]/.{f[x]->f,f'[x]->dfdx}/.(1+f)^(-1/2)->iwa
dwbdx1->D[(1+f[x] (((l1[x])^2)+((l2[x])^2)))^(1/2),x]/.
        {f[x]->f,f'[x]->dfdx1,l1[x]->l1,l2[x]->l2,l1'[x]->dl1dx1,l2'[x]->dl2dx1}/.
        {(1+f ((l1^2)+(l2^2)))^(-1/2)->iwb}
dwcdx1->D[(1+f[x] ((l2[x])^2))^(1/2),x]/.
        {f[x]->f,f'[x]->dfdx1,l1[x]->l1,l2[x]->l2,l1'[x]->dl1dx1,l2'[x]->dl2dx1}/.
        {(1+f (l2^2))^(-1/2)->iwc}

dtetraddx=FullSimplify[
D[Evaluate[
  tetrad/.{wa->wa[x,y,z],iwa->1/wa[x,y,z],
           wb->wb[x,y,z],iwb->1/wb[x,y,z],
           wc->wc[x,y,z],iwc->1/wc[x,y,z],
           f->f[x,y,z],
           l1->l1[x,y,z],l2->l2[x,y,z],l3->l3[x,y,z]}],x]//.
         
{Derivative[1,0,0][wa][x,y,z]->dwadx1,
 Derivative[1,0,0][wb][x,y,z]->dwbdx1,
 Derivative[1,0,0][wc][x,y,z]->dwcdx1,
 Derivative[1,0,0][f][x,y,z]->dfdx1,
 Derivative[1,0,0][l1][x,y,z]->dl1dx1,Derivative[1,0,0][l2][x,y,z]->dl2dx1,Derivative[1,0,0][l3][x,y,z]->dl3dx1,
 
 wa[x,y,z]->wa,
 wb[x,y,z]->wb,
 wc[x,y,z]->wc,
 f[x,y,z]->f,
 l1[x,y,z]->l1,l2[x,y,z]->l2,l3[x,y,z]->l3,
 
 wa^(-1)->iwa, wb^(-1)->iwb, wc^(-1)->iwc,
 wa^(-2)->iwasq, wb^(-2)->iwbsq, wc^(-2)->iwcsq}];

MatrixForm@dtetraddx[[f0]]
MatrixForm@dtetraddx[[f1]]
MatrixForm@dtetraddx[[f2]]
MatrixForm@dtetraddx[[f3]]

The derivatives of the tetrad should go to zero in the limit of M → 0.

MatrixForm[
Limit[Evaluate[dtetraddx//.{wa->(((1+f))^(1/2)),iwa->1/wa,iwasq->(iwa^2),
                            wb->(1+f ((l1^2)+(l2^2)))^(1/2),iwb->1/wb,iwbsq->(iwb^2),
                            wc->(1+f (l2^2))^(1/2),iwc->1/wc,
                            dwadx1->((dfdx1 iwa)/2),
                            dwbdx1->(1/2) iwb (f (2 dl1dx1 l1+2 dl2dx1 l2)+dfdx1 ((l1^2)+(l2^2))),
                            dwcdx1->(1/2) iwc (2 dl2dx1 f l2+dfdx1 (l2^2)),
                            f->2 M (r^3)/((r^4)+(a^2)(z^2)),
                            dfdx1->(((f^2) M qc x)/(2 qa (r^3)))}],{M->0}]]

Metric Derivatives

So now we have the tetrad, derivatives of the tetrad, the only missing ingredient is the derivatives of the metric.  
First find dr / dx and dr / dz:

(*dr/dx*)
FullSimplify[
FullSimplify[D[Evaluate[r/.{r->Sqrt[(((R^2)-(a^2)+Sqrt[((((R^2)-(a^2)))^2)+4(a^2)(z^2)])/2)]}/.{R->R[x]}],x]/.{R[x]->R,R'[x]->dRdx}]== 
FullSimplify[((dRdx R r  )/(Sqrt[((((R^2)-(a^2)))^2)+4 (a^2) (z^2)]))/.{r->Sqrt[(((R^2)-(a^2)+Sqrt[((((R^2)-(a^2)))^2)+4(a^2)(z^2)])/2)]}]]
Framed[drdx->FullSimplify[((dRdx R r  )/(Sqrt[((((R^2)-(a^2)))^2)+4 (a^2) (z^2)]))/.dRdx->x/R/.{(Sqrt[((((R^2)-(a^2)))^2)+4(a^2)(z^2)])^(-1)->(2 (r^2)+(a^2)-(R^2))^(-1)}]]

(*dr/dz*)
FullSimplify@PowerExpand[
FullSimplify[D[Evaluate[r/.{r->Sqrt[(((R^2)-(a^2)+Sqrt[((((R^2)-(a^2)))^2)+4(a^2)(z^2)])/2)]}/.{R->R[z]}],z]/.{R[z]->R,R'[z]->dRdz}]==
((dRdz R+(( 2dRdz (R^3)+(a^2) (-2 dRdz R+4 z))/(2(Sqrt[((((a^2)-(R^2)))^2)+4 (a^2) (z^2)]))))/(2 r))/.{r->Sqrt[(((R^2)-(a^2)+Sqrt[((((R^2)-(a^2)))^2)+4(a^2)(z^2)])/2)]}]
Framed[drdz->FullSimplify[FullSimplify[((dRdz R+(( 2dRdz (R^3)+(a^2) (-2 dRdz R+4 z))/(2(Sqrt[((((a^2)-(R^2)))^2)+4 (a^2) (z^2)]))))/(2 r))/.dRdz->z/R]/.
{(Sqrt[((((a^2)-(R^2)))^2)+4 (a^2) (z^2)])^(-1)->(2(r^2)-(R^2)+(a^2))^(-1)}]]

Now find df /dx and df / dz....(assuming m=1 here, as implemented in Kokkos version of Athena (for now) )

Framed[dfdx->FullSimplify[
FullSimplify[D[Evaluate[f/.{f->2(r^3)/((r^4)+(a^2)(z^2))}/.{r->r[x]}],x]/.{r[x]->r,r'[x]->drdx}] /. 
{(r^4)+(a^2)(z^2)->2(r^3)/f}]/.{-3(a^2)(z^2)->-qc-(r^4)}/.{drdx->((r x)/((a^2)+2 (r^2)-(R^2)))}/.{((a^2)+2 (r^2)-(R^2))->qa}]

dfdz->
FullSimplify[FullSimplify[FullSimplify[
D[Evaluate[f/.{f->2(r^3)/((r^4)+(a^2)(z^2))}/.{r->r[z]}],z]/.{r[z]->r,r'[z]->drdz}]/.
{-3(a^2)(z^2)->-qc-(r^4)}/.{(r^4)+(a^2)(z^2)->2(r^3)/f}] /.{drdz->((((a^2)+(r^2)) z)/(r ((a^2)+2 (r^2)-(R^2))))}/.((a^2)+(r^2))->qb];
FullSimplify[(((f^2) (-2 (a^2) r+((qb qc)/((a^2) r+2 (r^3)-r (R^2)))) z)/(2 (r^4)))==(((f^2) (-2 (a^2) (r^2)+((qb qc)/qa)) z)/(2 (r^5)))/.{qa->((a^2)+2 (r^2)-(R^2))}];
Framed[dfdz->(((f^2) (-2 (a^2) (r^2)+((qb qc)/qa)) z)/(2 (r^5)))]

Now find d l1/(dx^i) ,  d l2/(dx^i) ,  d l3/(dx^i)

FullSimplify[FullSimplify[FullSimplify[
D[(r x + a y)/((r^2) + (a^2))/.{r->r[x]},x]/.{r[x]->r,r'[x]->drdx}]/.
{drdx->((r x)/((a^2)+2 (r^2)-(R^2)))}]];
FullSimplify[x r((a^2)x-2a r y-(r^2)x)/((qb^2)qa)+r/(qb)==
((r (((a^2)+(r^2)) ((a^2)+2 (r^2)-(R^2)+(x^2))-2 r x (r x+a y)))/(((((a^2)+(r^2)))^2) ((a^2)+2 (r^2)-(R^2))))/.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl1dx->x r((a^2)x-2a r y-(r^2)x)/((qb^2)qa)+r/(qb)]

FullSimplify[FullSimplify[FullSimplify[
D[(r x + a y)/((r^2) + (a^2))/.{r->r[y]},y]/.{r[y]->r,r'[y]->drdy}]/.
{drdy->((r y)/((a^2)+2 (r^2)-(R^2)))}]];
FullSimplify[y r ((a^2)x-2 a r y-(r^2)x)/((qb^2)qa) + a/(qb)==
((-2 (r^2) y (r x+a y)+((a^2)+(r^2)) ((a^3)+2 a (r^2)-a (R^2)+r x y))/(((((a^2)+(r^2)))^2) ((a^2)+2 (r^2)-(R^2))))/.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl1dy->y r ((a^2)x-2 a r y-(r^2)x)/((qb^2)qa) + a/(qb)]

FullSimplify[FullSimplify[FullSimplify[
D[(r x + a y)/((r^2) + (a^2))/.{r->r[z]},z]/.{r[z]->r,r'[z]->drdz}]/.
{drdz->((((a^2)+(r^2)) z)/(r ((a^2)+2 (r^2)-(R^2))))}]];
FullSimplify[z/r ( (a^2) x - 2 a r y - (r^2) x )/(qb qa)==
((((a^2) x-(r^2) x-2 a r y) z)/(r ((a^2)+(r^2)) ((a^2)+2 (r^2)-(R^2)))) /.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl1dz->z/r ( (a^2) x - 2 a r y - (r^2) x )/(qb qa)]

FullSimplify[FullSimplify[FullSimplify[
D[(r y - a x)/((r^2) + (a^2))/.{r->r[x]},x]/.{r[x]->r,r'[x]->drdx}]/.
{drdx->((r x)/((a^2)+2 (r^2)-(R^2)))}]];
FullSimplify[x r ((a^2)y + 2 a r x - (r^2)y)/((qb^2)qa)-a/(qb)==
((-a ((a^2)+r (r-((2 r (x^2))/((a^2)+2 (r^2)-(R^2)))))+(((a-r) r (a+r) x y)/((a^2)+2 (r^2)-(R^2))))/((((a^2)+(r^2)))^2))/.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl2dx->x r ((a^2)y + 2 a r x - (r^2)y)/((qb^2)qa)-a/(qb)]

FullSimplify[FullSimplify[FullSimplify[
D[(r y - a x)/((r^2) + (a^2))/.{r->r[y]},y]/.{r[y]->r,r'[y]->drdy}]/.
{drdy->((r y)/((a^2)+2 (r^2)-(R^2)))}]];
FullSimplify[y r ((a^2)y + 2 a r x-(r^2)y)/((qb^2)qa) + r/(qb)==
((r (((a^2)+(r^2)) ((a^2)+2 (r^2)-(R^2))+2 a r x y+(a-r) (a+r) (y^2)))/(((((a^2)+(r^2)))^2) ((a^2)+2 (r^2)-(R^2))))/.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl2dy->y r ((a^2)y + 2 a r x-(r^2)y)/((qb^2)qa) + r/(qb)]

FullSimplify[FullSimplify[FullSimplify[
D[(r y - a x)/((r^2) + (a^2))/.{r->r[z]},z]/.{r[z]->r,r'[z]->drdz}]/.
{drdz->((((a^2)+(r^2)) z)/(r ((a^2)+2 (r^2)-(R^2))))}]];
FullSimplify[z/r ( (a^2) y + 2 a r x - (r^2) y )/(qb qa)==
(((2 a r x+(a^2) y-(r^2) y) z)/(r ((a^2)+(r^2)) ((a^2)+2 (r^2)-(R^2))))/.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl2dz->z/r ( (a^2) y + 2 a r x - (r^2) y )/(qb qa)]

FullSimplify[FullSimplify[FullSimplify[D[z/r/.{r->r[x]},x]/.{r[x]->r,r'[x]->drdx}]/.
{drdx->((r x)/((a^2)+2 (r^2)-(R^2)))}]];
FullSimplify[-x z/(r*qa)==-((x z)/((a^2) r+2 (r^3)-r (R^2)))/.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl3dx->-x z/(r*qa)]

FullSimplify[FullSimplify[FullSimplify[D[z/r/.{r->r[y]},y]/.{r[y]->r,r'[y]->drdy}]/.
{drdy->((r y)/((a^2)+2 (r^2)-(R^2)))}]];
FullSimplify[-y*z/(r*qa)==-((y z)/((a^2) r+2 (r^3)-r (R^2)))/.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl3dy->-y*z/(r*qa)]

FullSimplify[FullSimplify[FullSimplify[D[z/r/.{r->r[z]},z]/.{r[z]->r,r'[z]->drdz}]/.
{drdz->((((a^2)+(r^2)) z)/(r ((a^2)+2 (r^2)-(R^2))))}]];
FullSimplify[- (z^2)/(r^3) qb/qa +1/r ==(1/r)-((((a^2)+(r^2)) (z^2))/((r^3) ((a^2)+2 (r^2)-(R^2))))/.{qb->(r^2)+(a^2),qa->2(r^2)-(R^2)+(a^2)}]
Framed[dl3dz->- (z^2)/(r^3) qb/qa +1/r]

Now we have everything necessary to compute metric derivatives, e.g.,

((∂(Subscript[g, 00]))/(∂x))= ((∂f)/(∂x)) (Subscript[l, 0]) (Subscript[l, 0])    +    f ((∂l0)/(∂x)) (Subscript[l, 0])    +    f (Subscript[l, 0])  ((∂l0)/(∂x)).