I am attempting to create a Mathematica solver for the Einstein Field Equations (EFE) what will work for SI and Natural Units and any sign convention. It works for Natural Units, but when I try working with SI, I get results that don't agree with the literature. Specifically, for the metric tensor:$$g_{\mu\nu}=\begin{bmatrix}c^2 & 0 & 0 & 0\\0 & -\frac{a[t]^2}{1+k r^2} & 0 & 0\\0 & 0 & -r^2 a[t]^2 & 0\\0 & 0 & 0 & -r^2 a[t]^2 Sin[\theta]^2 \end{bmatrix}$$the solver yields this for the first solution:
$$\frac{c^2k}{a[t]^2}+\frac{a'[t]^2}{a[t]^2}=\frac{c^2\Lambda}{3}+\frac{8\pi G}{3}c^4\rho$$As you can see, it has a factor of $c^2$ in the density term that doesn't appear in the rest of the literature. I've traced the problem back to Covariant Stress Energy Momentum tensor for a perfect fluid:$$T_{\mu\nu}=\sum_{\alpha=1}^n g_{\mu\nu}T^{\mu}_{\nu}$$The time-time component of $T_{\mu\nu}$ is $T_{00}=c^4\rho$. This, in turn, can be traced back to the Trace of the Stress Energy Momentum Tensor:$$T^{\mu}_{\nu}=\begin{bmatrix}c^2\rho & 0 & 0 & 0\\0 & -p & 0 & 0\\0 & 0 & -p & 0\\0 & 0 & 0 & -p \end{bmatrix}$$The Stress Energy Momentum Tensor of a perfect fluid from [this article][1] is:$$T^{\mu\nu}=\left(\rho+\frac{p}{c^2}\right)u^\mu u^\nu-p g^{\mu\nu}$$ with $u=\{1,0,0,0\}$:$$T^{\mu\nu}=\begin{bmatrix}\rho & 0 & 0 & 0\\0 & \frac{p-k p r^2}{a[t]^2} & 0 & 0\\0 & 0 & \frac{p}{r^2 a[t]^2} & 0\\0 & 0 & 0 & \frac{p Csc[\theta]^2}{r^2 a[t]^2} \end{bmatrix}$$So the $c^2$ doesn't appear to come from the Stress-Energy Tensor. The only other option I see is that the term comes directly from the metric, $g_{00}=c^2$.
Could someone please tell me where the mistake is? This extra $c^4$ term in the first EFE solution appears to come directly from the metric tensor. What is the correct value for $T_{00}$ in SI units?


  [1]: https://en.wikipedia.org/wiki/Perfect_fluid