Hello Lubos:
This notebook sets out to detail the unified field proposal.  In effect, it serves as a "mathematical grammar checker."

The Lagrange density

These are the units of the components of the Lagrange density.

units = {γ -> 1, 
          J -> ((m^(1/2))/((L^(3/2)) t)), V -> (L^3), 
          q -> (((m^(1/2)) (L^(3/2)))/t), U -> (L/t), 
          (A^μ) -> ((m^(1/2))/(L^(1/2))), 
          (A^μ)ν -> ((m^(1/2))/(t (L^(1/2)))), (Sqrt[G]) -> 
            ((L^(3/2))/((m^(1/2)) t)), G -> ((L^3)/(m (t^2))), 
          c -> (L/t), h -> ((m (L^2))/t)};

Evaluate the units of the three parts of the Lagrange density: the kinetic energy, the mass and electric charge in motion, and the asymmetric field strength tensor.  If you are unfamiliar with Mathematica, "/." indicates rules for substitution, so /. units means the unit rules will be used to substitute into the preceding expression.  It is a way to check that all the parts of the Lagrange density have units of (m/(L^3)).

((J (A^μ))/c) /. units

((A^μ)(ν/(c^2)) /. units) (A^μ)ν /. units

Define the potential of the Lagrange density.

A = {ϕ[t, x, y, z], Ax[t, x, y, z], Ay[t, x, y, z], 
          Az[t, x, y, z]};

The following functions are required to check if the symmetric plus antisymmetric field strength tensor is identical with the asymmetric tensor.  One function flips the signs for covariant vectors. Four functions are required to take the contravariant and covariant derivative of a 4-vector with indices in uv or vu order.  A function is needed to contract two rank 2 tensors:

covariant[
        A_] := {A[[1]], -A[[2]], -A[[3]], -A[[4]]}

contraD[A_]:=(({{D[A[[1]], t], D[A[[2]], t], D[A[[3]], 
                t], D[A[[4]], t]}, {-c D[A[[1]], x], -c D[A[[2]], 
                  x], -c D[A[[3]], x], -c D[A[[4]], 
                  x]}, {-c D[A[[1]], y], -c D[A[[2]], 
                  y], -c D[A[[3]], y], -c D[A[[4]], 
                  y]}, {-c D[A[[1]], z], -c D[A[[2]], 
                  z], -c D[A[[3]], z], -c D[A[[4]], 
                  z]}}))

coD[A_]:=(({{D[A[[1]], t], D[A[[2]], t], D[A[[3]], 
                t], D[A[[4]], t]}, {c D[A[[1]], x], c D[A[[2]], x], c D[
                  A[[3]], x], c D[A[[4]], x]}, {c D[A[[1]], y], c D[A[[2]], y], c D[
                  A[[3]], y], c D[A[[4]], y]}, {c D[A[[1]], z], c D[A[[2]], z], c D[
                  A[[3]], z], c D[A[[4]], z]}}))

contraDvu[A_]:=(({{D[A[[1]], 
                t], -c D[A[[1]], x], -c D[A[[1]], 
                  y], -c D[A[[1]], z]}, {D[A[[2]], 
                t], -c D[A[[2]], x], -c D[A[[2]], 
                  y], -c D[A[[2]], z]}, {D[A[[3]], 
                t], -c D[A[[3]], x], -c D[A[[3]], 
                  y], -c D[A[[3]], z]}, {D[A[[4]], 
                t], -c D[A[[4]], x], -c D[A[[4]], 
                  y], -c D[A[[4]], z]}}))

coDvu[A_]:=(({{D[A[[1]], 
                t], c D[A[[1]], x], c D[A[[1]], 
                  y], c D[A[[1]], z]}, {D[A[[2]], 
                t], c D[A[[2]], x], c D[A[[2]], 
                  y], c D[A[[2]], z]}, {D[A[[3]], 
                t], c D[A[[3]], x], c D[A[[3]], 
                  y], c D[A[[3]], z]}, {D[A[[4]], 
                t], c D[A[[4]], x], c D[A[[4]], 
                  y], c D[A[[4]], z]}}))

contractMM[A_, B_] := 
      Sum[A[[i, j]] B[[i, j]], {i, 1, 4}, {j, 1, 4}]

Test the assertion that the contraction of the symmetric and antisymmetric tensors is equal to the contraction of the asymmetric tensor (will result in 0 if true):

Simplify[
      contractMM[contraD[A] + contraDvu[A], coD[covariant[A]] + coDvu[covariant[A]]]/(8 (c^2)) + 
        contractMM[contraD[A] - contraDvu[A], coD[covariant[A]] - coDvu[covariant[A]]]/(8 (c^2)) - 
        contractMM[contraD[A], coD[covariant[A]]]/(2 (c^2))]

Field equations

Define the electric and mass current densities:

(Subscript[J, q]) = {(Subscript[ρ, q]), (Subscript[J, qx]), (Subscript[J, qy]), (Subscript[J, qz])};

(Subscript[J, m]) = {(Subscript[ρ, m]), (Subscript[J, mx]), (Subscript[J, my]), (Subscript[J, mz])};

This package makes partial derivatives look like partial derivatives (it is not part of the standard set of packages).

<< FormatPD.m

Define the GEM Lagrange density:

Expand[(LGEM = -covariant[(Subscript[J, q]) - (Subscript[J, m])] . A/c - 
              Expand[contractMM[contraD[A], coD[covariant[A]]]/(2 (c^2))]) /. {ϕ[t, x, y, z] -> ϕ, Ax[t, x, y, z] -> Ax, 
          Ay[t, x, y, z] -> Ay, Az[t, x, y, z] -> Az}]

Define functions to apply the Euler-Lagrange equations to a Lagrange density.  The function potentialD takes the derivative of a Lagrange density with respect to the potential.  The function fieldD takes a second derivative with respect to the field of a Lagrange density.

potentialD[L_] := 
        Simplify[{D[L, ϕ[t, x, y, z]], D[L, Ax[t, x, y, z]], 
            D[L, Ay[t, x, y, z]], D[L, Az[t, x, y, z]]}];

fieldD[L_]:=Simplify[{D[D[L,Derivative[1, 0, 0, 0][ϕ][t, x, y, z]],t]+D[D[L,Derivative[0, 1, 0, 0][ϕ][t, x, y, z]],x]+D[D[L,Derivative[0, 0, 1, 0][ϕ][t, x, y, z]],y]+D[D[L,Derivative[0, 0, 0, 1][ϕ][t, x, y, z]],z],D[D[L,Derivative[1, 0, 0, 0][Ax][t, x, y, z]],t]+D[D[L,Derivative[0, 1, 0, 0][Ax][t, x, y, z]],x]+D[D[L,Derivative[0, 0, 1, 0][Ax][t, x, y, z]],y]+D[D[L,Derivative[0, 0, 0, 1][Ax][t, x, y, z]],z],D[D[L,Derivative[1, 0, 0, 0][Ay][t, x, y, z]],t]+D[D[L,Derivative[0, 1, 0, 0][Ay][t, x, y, z]],x]+D[D[L,Derivative[0, 0, 1, 0][Ay][t, x, y, z]],y]+D[D[L,Derivative[0, 0, 0, 1][Ay][t, x, y, z]],z],D[D[L,Derivative[1, 0, 0, 0][Az][t, x, y, z]],t]+D[D[L,Derivative[0, 1, 0, 0][Az][t, x, y, z]],x]+D[D[L,Derivative[0, 0, 1, 0][Az][t, x, y, z]],y]+D[D[L,Derivative[0, 0, 0, 1][Az][t, x, y, z]],z]}]

Apply to the GEM Lagrange density:

Simplify[Expand[c potentialD[LGEM]]] == 
      Expand[c fieldD[LGEM]]

Rewrite as a 4D wave equation.

(Subscript[J, q]) - (Subscript[J, m]) == (□^2) (A^μ);

Isolate Newton's field equation for gravity for the first 4D wave equation for the situation where q approaches zero and the second time derivative of phi is zero.

(c potentialD[LGEM][[1]] /. (Subscript[ρ, q]) -> 0)==(Expand[c fieldD[LGEM][[1]]]/. Derivative[2, 0, 0, 0][ϕ][t, x, y, z]->0)

Isolate Gauss' law for the first 4D wave equation where m approaches zero.

(c potentialD[LGEM][[1]] /. (Subscript[ρ, m]) -> 0) == 
      Expand[c fieldD[LGEM][[1]]]

Classical fields

The long name "EField" for E must be used since E means 2.718... to Mathematica.  Define the five classical fields that constitute the asymmetric tensor A^u,v

Efield = {-D[Ax[t, x, y, z], t] - 
          D[ϕ[t, x, y, z], x], -D[Ay[t, x, y, z], t] - 
          D[ϕ[t, x, y, z], y], -D[Az[t, x, y, z], t] - 
          D[ϕ[t, x, y, z], z]}
e = {D[Ax[t, x, y, z], t] - D[ϕ[t, x, y, z], x], 
        D[Ay[t, x, y, z], t] - D[ϕ[t, x, y, z], y], 
        D[Az[t, x, y, z], t] - 
          D[ϕ[t, x, y, z], z]}
B = 
      Curl[{Ax[t, x, y, z], Ay[t, x, y, z], 
          Az[t, x, y, z]}]
b = {-D[Ay[t, x, y, z], z] - 
          D[Az[t, x, y, z], y], -D[Ax[t, x, y, z], z] - 
          D[Az[t, x, y, z], x], -D[Ax[t, x, y, z], y] - 
          D[Ay[t, x, y, z], x]}
g = {D[ϕ[t, x, y, z], 
          t], -D[Ax[t, x, y, z], x], -D[Ay[t, x, y, z], y], -D[
            Az[t, x, y, z], z]}

These functions are needed to determine the components of the asymmetric and symmetric field strength tensors.

contraDvu[A_]:=(({{D[A[[1]], 
                t], -c D[A[[1]], x], -c D[A[[1]], 
                  y], -c D[A[[1]], z]}, {D[A[[2]], 
                t], -c D[A[[2]], x], -c D[A[[2]], 
                  y], -c D[A[[2]], z]}, {D[A[[3]], 
                t], -c D[A[[3]], x], -c D[A[[3]], 
                  y], -c D[A[[3]], z]}, {D[A[[4]], 
                t], -c D[A[[4]], x], -c D[A[[4]], 
                  y], -c D[A[[4]], z]}}))
symmetricD[A_] := 
      contraD[A] + contraDvu[A]
antisymmetricD[
        A_] := contraD[A] - contraDvu[A]

Write out the antisymmetric (E + B), symmetric (e + b + g), and a symmetric tensors (all) in terms of the individual components.

MatrixForm[antisymmetricD[A]]

MatrixForm[symmetricD[A]]

\(\(\(
\)\(Simplify[
      MatrixForm[\((symmetricD[A] + antisymmetricD[A])\)/2]]\)\)\)

To do div, grad, curl and all that requires the vector analysis package:

<< Calculus`VectorAnalysis`

SetCoordinates[Cartesian[x, y, z]];

GEM version of a unified Gauss' law:

Simplify[(c/2) (Div[Efield] + Div[e]) + 
          D[g[[1]], t]/c] == ((Subscript[J, q]) - (Subscript[J, m]))[[1]]

These sets of substitution rules are required to set DivE = 0 and Div e = 0:

noEfield = {D[Ax[t, x, y, z], t] -> D[ϕ[t, x, y, z], x], 
          D[Ay[t, x, y, z], t] -> D[ϕ[t, x, y, z], y], 
          D[Az[t, x, y, z], t] -> D[ϕ[t, x, y, z], z]};

noe = {D[Ax[t, x, y, z], t] -> -D[ϕ[t, x, y, z], x], 
          D[Ay[t, x, y, z], t] -> -D[ϕ[t, x, y, z], y], 
          D[Az[t, x, y, z], 
              t] -> -D[ϕ[t, x, y, z], z]};

nogt = {D[ϕ[t, x, y, z], t] -> 0};

If there is no divergence of the E field, no dynamic g, and no electric charge density, Newton's field equations for gravity results.

Simplify[(c/2) (Div[Efield /. noEfield] + Div[e /. noEfield]) + 
            D[g[[1]] /. nogt]/c] == ((Subscript[J, q]) - 
              (Subscript[J, m]))[[1]] /. (Subscript[ρ, q]) -> 0

If there is no divergence of the symmetric e filed and m is zero, Gauss' law results.

Simplify[(c/2) (Div[Efield /. noe] + Div[e /. noe]) + 
            D[g[[1]], t]/c] == ((Subscript[J, q]) - 
              (Subscript[J, m]))[[1]] /. (Subscript[ρ, m]) -> 0

The homogeneous Maxwell equations

Div[B]

Curl[Efield] + D[B, t]

A physically relevant 4D solution

A function to test for solutions to the vacuum field equations:

test[potential_] := 
      Simplify[{D[potential[[1]], {t, 2}] - 
            D[potential[[1]], {x, 2}] - D[potential[[1]], {y, 2}] - 
            D[potential[[1]], {z, 2}], 
          D[potential[[2]], {t, 2}] - D[potential[[2]], {x, 2}] - 
            D[potential[[2]], {y, 2}] - D[potential[[2]], {z, 2}], 
          D[potential[[3]], {t, 2}] - D[potential[[3]], {x, 2}] - 
            D[potential[[3]], {y, 2}] - D[potential[[3]], {z, 2}], 
          D[potential[[4]], {t, 2}] - D[potential[[4]], {x, 2}] - 
            D[potential[[4]], {y, 2}] - 
            D[potential[[4]], {z, 2}]}]

Test the 1/R potential. The units of a potential are ((Sqrt[m])/(Sqrt[L])).

test[{((Sqrt[(h/c)])/(Sqrt[(x^2) + (y^2) + (z^2)])), 0, 0, 0}]

(1/L) (Sqrt[(h/c)]) /. units

Test a potential that is an inverse distance squared.

test[{(((Sqrt[G]) h/(c^2))/((x^2) + (y^2) + (z^2) - (t^2))), 0, 0, 
        0}]

(((Sqrt[G]) h)/((c^2) (L^2))) /. units

The units imply this potential would involve relativistic quantum gravity.

Test a potential that is a normalized linear perturbation.

a = {(((c/(Sqrt[G])))/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))), 0, 0, 0};

test[a]

(c/(Sqrt[G])) /. units

These units imply relativistic gravity, but not quantum mechanics.  Variations on this potential are the focus of the following work.

Look at the derivatives to first order in k.

Expand[D[a, t]] /. (k^2) -> 0

For a small oscillation, the denominator will be approximately 2  (σ^2).  This substitution list will make it so (and for other cases to come).

sublist = {(k^2) -> 
            0, -(((1 + ((k t)/(σ^2))))^2) + (((1 + ((k x)/(σ^2))))^2) + (((1 + ((k y)/(σ^2))))^2) + (((1 + ((k z)/(σ^2))))^2) -> 
            2, -(((1 + ((k t)/(σ^2))))^2) + (((1 - ((k x)/(σ^2))))^2) + (((1 - ((k y)/(σ^2))))^2) + (((1 - ((k z)/(σ^2))))^2) -> 
            2, -(((1 - ((k t)/(σ^2))))^2) + (((1 + ((k x)/(σ^2))))^2) + (((1 - ((k y)/(σ^2))))^2) + (((1 - ((k z)/(σ^2))))^2) -> 
            2, -(((1 - ((k t)/(σ^2))))^2) + (((1 - ((k x)/(σ^2))))^2) + (((1 + ((k y)/(σ^2))))^2) + (((1 - ((k z)/(σ^2))))^2) -> 
            2, 
-(((1 - ((k t)/(σ^2))))^2) + (((1 - ((k x)/(σ^2))))^2) + (((1 - ((k y)/(σ^2))))^2) + (((1 + ((k z)/(σ^2))))^2) -> 
            2, 
-(((1 - ((k t)/(σ^2))))^2) + (((1 + ((k x)/(σ^2))))^2) + (((1 + ((k y)/(σ^2))))^2) + (((1 + ((k z)/(σ^2))))^2) -> 
            2, 
-(((1 - ((k t)/(σ^2))))^2) + (((1 - ((k x)/(σ^2))))^2) + (((1 - ((k y)/(σ^2))))^2) + (((1 - ((k z)/(σ^2))))^2) -> 2};

Expand[D[a, t]] /. sublist

Write a potential that solves the field equations, but only has derivatives along the diagonal.

diagonalSHO = {(1/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))) + 
            (1/((((1 - (((k x)/(σ^2)))))^2) + (((1 - (((k y)/(σ^2)))))^2) + (((1 - (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))), 
          (1/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))) + 
            (1/((((1 + (((k x)/(σ^2)))))^2) + (((1 - (((k y)/(σ^2)))))^2) + (((1 - (((k z)/(σ^2)))))^2) - (((1 - (((k t)/(σ^2)))))^2))), 
          (1/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))) + 
            (1/((((1 - (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 - (((k z)/(σ^2)))))^2) - (((1 - (((k t)/(σ^2)))))^2))), 
          (1/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))) + 
            (1/((((1 - (((k x)/(σ^2)))))^2) + (((1 - (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 - (((k t)/(σ^2)))))^2)))};

test[diagonalSHO]

diagonalSHO2 = {((-1)/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 - (((k t)/(σ^2)))))^2))) + ((-1)/((((1 - (((k x)/(σ^2)))))^2) + (((1 - (((k y)/(σ^2)))))^2) + (((1 - (((k z)/(σ^2)))))^2) - (((1 - (((k t)/(σ^2)))))^2))), 
          (1/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))) + 
            (1/((((1 + (((k x)/(σ^2)))))^2) + (((1 - (((k y)/(σ^2)))))^2) + (((1 - (((k z)/(σ^2)))))^2) - (((1 - (((k t)/(σ^2)))))^2))), 
          (1/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))) + 
            (1/((((1 - (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 - (((k z)/(σ^2)))))^2) - (((1 - (((k t)/(σ^2)))))^2))), 
          (1/((((1 + (((k x)/(σ^2)))))^2) + (((1 + (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 + (((k t)/(σ^2)))))^2))) + 
            (1/((((1 - (((k x)/(σ^2)))))^2) + (((1 - (((k y)/(σ^2)))))^2) + (((1 + (((k z)/(σ^2)))))^2) - (((1 - (((k t)/(σ^2)))))^2)))};

test[diagonalSHO2]

diagonalSHO2 /. k -> 0

Take the derivative of the potential.

MatrixForm[Expand[contraD[diagonalSHO]] /. sublist]

MatrixForm[Expand[contraD[diagonalSHO2]] /. sublist]

Use this derivative of a potential in a Lorentz-like force law

contractVM[V_, M_] := {Sum[V[[i]] M[[1, i]], {i, 1, 4}], 
        Sum[V[[i]] M[[2, i]], {i, 1, 4}], 
        Sum[V[[i]] M[[3, i]], {i, 1, 4}], 
        Sum[V[[i]] M[[4, i]], {i, 1, 4}]}

v = {(Subscript[U, 0])[τ], (Subscript[U, 1])[τ], (Subscript[U, 2])[τ], 
          (Subscript[U, 3])[τ]};

Simplify[-m c contractVM[covariant[v], 
          Expand[contraD[diagonalSHO]] /. sublist]]

Replace (σ^2)with -(τ^2).  From here on, assume spherical symmetry, so (Subscript[U, R]) does the work of (Subscript[U, 1]),(Subscript[U, 2]), and (Subscript[U, 3]).

-1*{-((c k m (Subscript[U, 0])[τ])/(σ^2)), ((c k m (Subscript[U, 1])[τ])/(σ^2))} /. σ -> c τ

Apply the chain rule to the change in momentum, and assume that dm/dtau is zero.

{D[m[τ] (Subscript[U, 0])[τ, n], τ], 
        D[m[τ] (Subscript[U, R])[τ, n], τ]}/.m'[τ]->0

Solve the first order differential equations.

DSolve[{((k m (Subscript[U, 0])[τ])/(c (τ^2))) == 
          m D[(Subscript[U, 0])[τ], τ], -((k m (Subscript[U, R])[τ])/(c (τ^2))) == m D[(Subscript[U, R])[τ], τ]}, {(Subscript[U, 0])[τ], 
        (Subscript[U, R])[τ]}, τ]

Contract the relativistic velocity 4-vector solution.

{(E^(-(k/(c τ)))) C[
            2], (E^((k/(c τ)))) C[
            1]} . {(E^(-(k/(c τ)))) C[
            2], -(E^((k/(c τ)))) C[1]}

In flat spacetime, k->0

-(E^(((2 k)/(c τ)))) C[1]^2 + (E^(-((2 k)/(c τ)))) C[2]^2 == (c^2) /. 
      k -> 0

This can only be the case if the constants of integration are velocities in flat spacetime.

flatSpacetimeConstraint = { C[1] -> (dR/dtau), 
          C[2] -> c (dt/dtau)};

Apply this constraint to the contracted relativistic velocity 4-vector solution.

-(E^(((2 k)/(c τ)))) C[1]^2 + (E^(-((2 k)/(c τ)))) C[2]^2 == (c^2) /. 
      flatSpacetimeConstraint

Multiply both sides by ((dtau^2)/(c^2)).

Expand[((((c^2) (dt^2) (E^(-((2 k)/(c τ)))))/(dtau^2)) - (((dR^2) (E^(((2 k)/(c τ)))))/(dtau^2))) ((dtau^2)/(c^2))] == (c^2) ((dtau^2)/(c^2))

To apply to a weak gravitational system, substitute τ -> (R/c) and k -> (GM/(c^2)).

((dt^2) (E^(-((2 k)/(c τ)))) - (((dR^2) (E^(((2 k)/(c τ)))))/(c^2)) /. {τ -> 
              (R/c),  k -> ((G M)/(c^2))}) == (dtau^2)

This is the exponential metric which is consistent with experimental tests of general relativity.

Added Note

Consider a static electrically neutral system in a vacuum.  The Gauss-like law for this proposal is this:

▽e  =  ▽ ((A^(m, ν)) + 
              (A^(ν, μ)) - (((Subscript[Γ, ω]) )^i0)) = 0

Choose a reference frame such that the derivative of the potential happens to be zero.   Under these condistions, we want to calculate the divergence of the Christoffel symbol.

((▽(Subscript[Γ, ω]) )^i0) = ▽(Subscript[g, ω β]) ((g^(β i,  0)) + (g^(0 β,  i)) - 
            (g^(i 0,  β))) (A^ω)

The first term drops because the metric is static. The third term drops becuase the metric is diagonal. Beta must equal 0 for the metric to be non-zero.

g00 = 
        Exp[-((2 G M)/((c^2) (Sqrt[(x^2) + (y^2) + (z^2)])))];
(Subscript[g, 00]) = 
        Exp[((2 G M)/((c^2) (Sqrt[(x^2) + (y^2) + (z^2)])))];
g11 = g22 = g33 = 
            Exp[((2 G M)/((c^2) (Sqrt[(x^2) + (y^2) + (z^2)])))];

(Subscript[g, 00]) g00

Singular solutions

Simplify[-((c^2)/(Sqrt[G])) (D[
              (Subscript[g, 00]) D[g00, x], x] + D[(Subscript[g, 00]) D[g00, y], y] + 
            D[(Subscript[g, 00]) D[g00, z], z])]