Compute curvature of the Taub-NUT metric

taubNUT.m

rhosq = r^2 + nut^2
Deltasq = r^2 - 2 M r - nut^2

(* Covariant metric tensor (metric on vectors) *)
g = Simplify[{{-Deltasq/rhosq,0,0,-2 nut Cos[theta] Deltasq/rhosq},{0,rhosq/Deltasq,0,0},{0,0,rhosq,0},{-2 nut Cos[theta] Deltasq/rhosq,0,0,rhosq Sin[theta]^2 - 4 nut^2 Cos[theta]^2 Deltasq/rhosq}}]

(* Contravariant metric tensor (metric on covectors) *)
ginv = Simplify[Inverse[g]]

coords = {t,r,theta,phi}

(* First kind Christoffel symbols [i j, k] *)
chris1 = Simplify[Table[(-D[g[[i,j]],coords[[k]]]+D[g[[i,k]],coords[[j]]]+D[g[[j,k]],coords[[i]]])/2,{i,1,4},{j,1,4},{k,1,4}]]

(* Second kind Christoffel symbols {l \\ i k} *)
chris2 = Simplify[Table[Sum[chris1[[i,j,k]] ginv[[k,l]],{k,1,4}],{l,1,4},{i,1,4},{j,1,4}]]

(* Riemann tensor R^l_ijk *)
riemann = Simplify[Table[D[chris2[[l,i,k]],coords[[j]]]-D[chris2[[l,i,j]],coords[[k]]]+Sum[chris2[[m,i,k]] chris2[[l,m,j]],{m,1,4}]-Sum[chris2[[m,i,j]] chris2[[l,m,k]],{m,1,4}],{l,1,4},{i,1,4},{j,1,4},{k,1,4}]]

(* Ricci tensor *)
(* Check: this should be zero: *)
ricci = Simplify[Table[Sum[riemann[[l,i,j,l]],{l,1,4}],{i,1,4},{j,1,4}]]

(* Fully covariant metric tensor R_lijk *)
riemanncov = Simplify[Table[Sum[riemann[[m,i,j,k]] g[[l,m]],{m,1,4}],{l,1,4},{i,1,4},{j,1,4},{k,1,4}]]

(* Kretschmann curvature scalar *)
kretschmann = Simplify[Sum[riemanncov[[l1,i1,j1,k1]] riemanncov[[l2,i2,j2,k2]] ginv[[l1,l2]] ginv[[i1,i2]] ginv[[j1,j2]] ginv[[k1,k2]], {l1,1,4},{l2,1,4},{i1,1,4},{i2,1,4},{j1,1,4},{j2,1,4},{k1,1,4},{k2,1,4}]]
Series[kretschmann, {r,Infinity,6}]

See https://bsky.app/profile/gro-tsen.bsky.social/post/3l5eh53ye6z2f for context for this computation.