In this notebook, I will propose a unified field Lagrange density to unify gravity and EM.  Here is the list of topics to be covered:
    I:   The Lagrange density itself
    II:  Field equations
    III: Classical fields
    IV: Dynamic metrics
    V: New constant velocity solutions
    VI. The Stress-energy tensor
    VII:  Quantization
Everything will be done using standard techniques in four dimensions.

I: The GEM Lagrange density

What is a Lagrange density? It is a scalar function that describes all the mass and energy interactions per unit volume, bar none. The complete description of all energy as it happens to be distributed per unit volume is why a Lagrange density is important. Once the Lagrange density is set, everything else follows, including things like the field equations, energy/momentum, and whether the proposal can be quantized. The logic of mathematical physics is unbending.

The core idea is to start with the classical Maxwell lagrange density, and generize it enough so that the Lagrangian can also describe gravity. Here is the Lagrange density that will be studied:

ℒ = -(1/c) ((Subscript[J, q]^μ) - (Subscript[J, m]^μ)) 
          (Subscript[A, μ]) - (1/(2 
                 (c^2))) (∇^μ) (A^ν) (Subscript[∇, μ]) (Subscript[A, ν])
where
ℒ = The Lagrange density
(Subscript[J, q]^μ) = 
      electric charge contravariant 4 - current density 
(Subscript[J, m]^μ) = 
      mass charge contravariant 4 - 
        current density 
(Subscript[A, μ]) = gravity/EM covariant 4 - potential
(∇^μ) (A^ν) = 
      contravariant 4 - derivative of a contravariant 4 - 
        potential

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

units = {J -> ((m^(1/2))/((L^(3/2)) t)), V -> (L^3), 
          q -> (((m^(1/2)) (L^(3/2)))/t), U -> (L/t), 
          Amu -> ((m^(1/2))/(L^(1/2))), 
          Amunu -> ((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 Lagrange density: 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 Amu)/c) /. units

((Amunu/(c^2)) /. units) Amunu /. units

The charge coupling term (J.A) represents the energy of the electric and mass charges in motion in a potential per unit volume.  The field strength tensor term represents all the energy in a changing potential in flat or curved spacetime. The idea of treating mass technically as a charge is not new, but is also not widespread.  One basic question is how to get the units such that electric and mass charge are identical.  In the Gaussian system, the units for q can be found from Coulomb's law: F = ((q q')/(R^2)).  The units of force are ((m L)/(t^2)), so (q^2) = ((m (L^3))/((t^2))), or q = ((Sqrt[m (L^3)])/t).  Odd, but true.  Determine the units of (Sqrt[G]) m: (Sqrt[(((L^3))/(m(t^2)))]) m = ((Sqrt[m (L^3)])/t) so the units of q and (Sqrt[G]) m are identical. "(Sqrt[G]) m" will be referred to as "mass charge" so that the language parallels electric charge.

Half of the terms in this Lagrange density are part of the standard classical Lagrange density for EM: the electric current density coupling to the potential field (-(Subscript[J, q]^μ)  (Subscript[A, μ])) and the antisymmetric electromagnetic field strength tensor ((F^μ)ν=(∂^μ)(A^ν)- (∂^ν)(A^μ)) which gets contracted with itself.

One new term is a mass current density coupling to the same potential field (+(Subscript[J, m]^μ)  (Subscript[A, μ])).  One might wonder why (Subscript[J, q]^μ) is not redefined to include the mass current density.  The reason has to do with the sign difference between the current density coupling and the field strength tensor.  For EM, the same sign for both the EM current coupling and field strength tensor will lead to field equations where like electrical charges repel.  For the mass charge, the difference in the sign between the mass current coupling and the field strength tensor will lead to field equations where like mass charges attract, a key property of gravity.

The second additional term is the contraction of a symmetric field strength tensor ((∇^μ)(A^ν)+ (∇^ν)(A^μ)).  In order to transform like a tensor, it requires covariant derivatives ("∇") unlike the antisymmetric tensor's exterior derivative ("∂"). The covariant derivative is a way to say the 4-derivative depends on how a metric changes.  The exterior derivative is unaffected by a metric changing.  Why the difference?  A metric is a symmetric tensor.  Any changes in the metric will only show up in the symmetric tensor.

Let's think about the symmetry of the field strength tensors.  One type of symmetry is about the order of the indices.  If μ changes places with ν, the symmetric tensor ((∇^μ)(A^ν)+ (∇^ν)(A^μ)) will not be changed, but all the signs in the antisymmetric tensor ((∂^μ)(A^ν)- (∂^ν)(A^μ)) would have to be flipped. This may be why there is one charge for mass and two for electricity.  The asymmetric tensor  (∇^μ) (A^ν) does not have a simple relationship with (∇^ν) (A^μ).  It is also called reducible because it can be represented by the irreducible symmetric and antisymmetric field strength tensors.

Because there is no t, x, or angles in the Lagrangian, energy, linear, and angular momentum will be conserved. This is a technical advantage over the Hilbert action of general relativity where energy cannot be defined locally. EM is a gauge theory, as is general relativity. Here too there is a gauge, a choice one must make before doing a calculation. The choice of gauge arrises from the definition of a covariant derivative: (∇^μ)(A^ν)=\(∂\^μ A\^ν - \(Γ\_σ\^μν\) A\^σ\). For a manifold with a torsion-free, metric compatible connection, there is one metric that uniquely determines the connection, as is done in general relativity. One has the ability to choose a gauge which is about the relationship between changes in the potential and changes in the metric. One could choose to work in flat Euclidean spacetime, in which case all the change found in the covariant derivative will be due to the standard derivative term, (∂^μ) (A^ν). Or the potential could be chosen so it makes no contribution to the contravariant derivative, and everything depends on the change in metric, ((Subscript[Γ, σ])^μ)ν 
        (A^σ). There are an infinite number of gauge choices between these two extremes. The potential/connection diffeomorphism is the key unifying link between gravity and EM.

Even though I worked with asymmetric, symmetric, and antisymmetric tensors for a while, I still wondered about their relationships to each other because I have little day-to-day experience with indices. One way to see precisely what is going on is to start with an asymmetric tensor represented as a 4x4 real matrix filled with random integers:

B=(({{1, 1, 2, 3}, {5, 6, 9, 8}, {0, 7, 11, 12}, {13, 0, 16, 0}}));

Take the transpose of this matrix to get (B^ν)μ, and it should be obvious: (B^μ)ν does not have any simple relationship to (B^ν)μ.  The challenge is to rewrite this matrix as the sum of a symmetric and an antisymmetric matrix.  Try it!

Here is the simple idea: the symmetric matrix is the average values of the two values on opposite sides of the diagonal, and the antisymmetric matrix is the deviation from those average values:

MatrixForm[BAve  = ((B + Transpose[B])/2)]

MatrixForm[BDev  = ((B - Transpose[B])/2)]

B - (BAve  +  BDev)

Isn't that cool?  Taking a transpose of the first matrix is not going to change an average value.  Taking a transpose of the deviations from average will flip all the signs, but no magnitudes.  The average amount of change symmetric matrix takes the average values off the diagonal and the diagonal itself.  The symmetric matrix represents ten of the sixteen terms of the asymmetric tensor.  For the antisymmetric matrix, the diagonal is all zeroes, and represents the remaining six terms. The average value can be independent of the deviation from the average.

Now we can use new words to describe the symmetric field strength tensor and the antisymmetric tensor.  The average amount of change in the 4-potential is (∇^μ)(A^ν)+ (∇^ν)(A^μ), which includes changes that occur due to a changing metric.  How much the metric changes in involved in the gauge choice.  The deviation from that average amount of change is the tensor (∂^μ)(A^ν)- (∂^ν)(A^μ), which excludes any contribution by a changing metric.  The complete covariant 4-derivative of a contravariant 4-potential field strength tensor is the average amount of change symmetric tensor plus the deviation from the average amount of change antisymmetric tensor.

In summary, the Lagrange density under study has two parts: a coupling term for electric and mass charge in motion in a potential, and the average amount of change and deviation from the average amount of change in a 4-potential field strength tensors.

II: Field equations

What does one do with a Lagrange density?  It cannot be measured directly.  Instead, by taking certain derivatives, things which are physically observable are found.  In this section the field equations are derived.  Field equations are used to describe the motion and distribution of particles in a volume of spacetime.

How does one generate field equations?  Apply the Euler-Lagrange equation to a Lagrange density.  It is that formulaic.

It would be better to understand how the Euler-Lagrange equation works its magic.  Here is a sketch.  Consider what would happen if the Lagrange density were integrated over a volume and time in curved spacetime:

𝒮 = ∫ℒ (Sqrt[-det g]) 
          dV dt

The action is the integral of all the mass and energy interactions per unit volume over a volume and time.  If ℒ is varied, but the value of the integral did not change over arbitrary times, then both a symmetry of the action and a conserved quantity has been found.

If the Lagrange density depends only on the potential (A^ν) and the field strength tensor (Subscript[∇, μ]) (A^ν) because everything else is known up to a diffeomorphism, the variation of the action 𝒮 will always be zero if:

∂ (ℒ/∂) (A^ν)=(∇^μ)(((∂ ℒ)/(∂((∇^μ)(A^ν)))))

This is the Euler-Lagrange equation.

To get Mathematica to do this 19th century math requires that a few things get defined:

covariantvec[
        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]}}))
symmetricD[A_] := 
      contraD[A] + contraDvu[A]
antisymmetricD[
        A_] := contraD[A] - contraDvu[A]

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

Define the gravity/EM potential, and the current densities for electric and mass charges:

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

(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])};

Now we can define all the terms in the GEM Lagrangian:

Expand[(LGEM = -covariantvec[(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.  [Notation is needed so the output looks intelligible to people.]

<< Utilities`Notation`

Notation[\(∂f_\/∂t_\) ⟺ Derivative[1][f_][t_]]
Notation[\(∂\^n_ f_\/∂t_\^n_\) ⟺ Derivative[n_][f_][t_]]
Notation[ ∂ f_[t_, x_, y_, z_]/∂t_ ⟺ Derivative[1, 0, 0, 0][f_][t_, x_, y_, z_]]
Notation[ ∂ f_[t_, x_, y_, z_]/∂x_ ⟺ Derivative[0, 1, 0, 0][f_][t_, x_, y_, z_]]
Notation[ ∂ f_[t_, x_, y_, z_]/∂y_ ⟺ Derivative[0, 0, 1, 0][f_][t_, x_, y_, z_]]
Notation[ ∂ f_[t_, x_, y_, z_]/∂z_ ⟺ Derivative[0, 0, 0, 1][f_][t_, x_, y_, z_]]
Notation[\(\(\ \)\(∂\^n_\ 
                f_[t_, x_, y_, z_]\/∂t_\^n_\)\) ⟺ Derivative[n_, 0, 0, 0][f_][t_, x_, y_, z_]]
Notation[\(\(\ \)\(∂\^n_\ 
                    f_[t_, x_, y_, z_]\/∂x_\^n_\)\) ⟺ Derivative[0, n_, 0, 0][f_][t_, x_, y_, z_]]Notation[\(\(\ \)\(∂\^n_\ 
                    f_[t_, x_, y_, z_]\/∂y_\^n_\)\) ⟺ Derivative[0, 0, n_, 0][f_][t_, x_, y_, z_]]Notation[\(\(\ \)\(∂\^n_\ 
                    f_[t_, x_, y_, z_]\/∂z_\^n_\)\) ⟺ Derivative[0, 0, 0, n_][f_][t_, x_, y_, z_]];

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]]) /. {ϕ[t, x, y, 
            z] -> ϕ}

This is a 4D wave equation with two charges, one where like electric charges repel and another where like mass charges attract.

Isolate Gauss' static law for the first 4D wave equation where the mass charge density approaches zero.

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

Isolate Newton's field equation for gravity for the first 4D wave equation for the situation where the electric charge density goes to 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)

This the mass charge density, not mass density. This must be the case for the units to be consistent using the same gravity/EM potential. It is a fun exercise to apply these two equations to a proton point charge at 1cm.  The mass charge is something like 13 orders of magnitude smaller than the electric charge, and we only can measure the electric charge to ten significant digits.

"General Gauss' law" is my name for the covariant unified field law modeled on the one from EM:

-c potentialD[LGEM][[1]] == -Expand[
          c fieldD[LGEM][[1]]]

The covariant form is expected for Gauss' law, but is significant news for Newton's law of gravity.

III: 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^ν):

Efield = {-D[Ax[t, x, y, z], t] - 
          c D[ϕ[t, x, y, z], x], -D[Ay[t, x, y, z], t] - 
          c D[ϕ[t, x, y, z], y], -D[Az[t, x, y, z], t] - 
          c D[ϕ[t, x, y, z], z]}
e = {D[Ax[t, x, y, z], t] - c D[ϕ[t, x, y, z], x], 
        D[Ay[t, x, y, z], t] - c D[ϕ[t, x, y, z], y], 
        D[Az[t, x, y, z], t] - 
          c D[ϕ[t, x, y, z], z]}
B = 
      c Curl[{Ax[t, x, y, z], Ay[t, x, y, z], 
            Az[t, x, y, z]}]
b = -c {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], -c D[Ax[t, x, y, z], x], -c D[Ay[t, x, y, z], 
            y], -c D[Az[t, x, y, z], z]}

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]];

To get the signs correct, it is vital to note one is using a contravariant derivative ((∇^μ)) to contract with the field strength tensor ((Subscript[∇, ν]) (A^ν)). The operators contraDiv and contraCurl take this into account. Here is the GEM version of a unified Gauss' law:

symcurl[{a_, b_, c_}] := {D[b, z] + D[c, y], D[c, x] + D[a, z], 
        D[a, y] + D[b, x]}
grad3[{a_, b_, c_, d_}] := {D[b, x], D[c, y], D[d, z]}
contraGrad3[v_] := -grad3[v]

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

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.

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

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

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

Ampere's Law:

{((Subscript[J, q]) - (Subscript[J, m]))[[2]], ((Subscript[J, q]) - (Subscript[J, m]))[[3]], ((Subscript[J, q]) - 
              (Subscript[J, m]))[[4]]} == 
      Simplify[(1/2) (D[-Efield, t]/c + D[e, t]/c + Curl[B] + 
                symcurl[b]) + grad3[g]]

The homogeneous Maxwell equations

Simplify[Div[B]]

Simplify[Curl[Efield] + D[B, t]/c]

There are no gravitational analogs to the homogeneous Maxwell equations.

IV: Dynamics Metrics

All of the equations written so far have been manifestly covariant, except for the static laws. This means that their form will not change no matter what metric is used. The equations will not change their form no matter how the metric changes throughout the spacetime manifold. The metric may or may not satisfy the Einstein field equations: either way, the field equations remain the same. Although it is natural to presume all the field equations were written with a flat Minkowski metric, technically that does not have to be the case.

In this section, I will show two different roads to a dynamic metric which is consistent with experimental tests of weak field gravity. The first approach looks at Gauss' law when the derivatives are treated as covariant derivatives. The second road to the identical metric involves finding a solution to the wave equation, plugging that into a force law, and rearranging the result to look like a metric.

Consider a static, spherically symmetric, system in a vacuum.  The Gauss-like law for this proposal is this:

\(\(1\/2\) \((∇\(\(.\)\(E\)\) + ∇\(\(.\)\(e\)\))\) = \(∇\_μ\((∂\^μ A\^0 + ∂\^0 A\^μ - 
              2 \( Γ\_σ\^μ0\) A\^σ)\) = 
        0\)\)

Choose a potential such that the derivative of the potential happens to be zero. This is effectively a choice of gauge so that the dynamic metric contains all the information about the mass and electric charges in the system.   Under these conditions, calculate the divergence of the Christoffel symbol.

\(\(\(\(▽Γ\_σ\)\(\ \)\)\^i0\) 
        A\^σ = \(1\/2\) ▽\ \(g\_\(σ\ β\)\) \((g\^\(β\ i, 0\) + g\^\(0\ β, i\) - g\^\(i\ 0, β\))\) 
        A\^σ\)

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

Now we need to "guess" a metric that will solve this equation. The metric must reduce to the Minkowski metric for a small mass. It must solve a Poisson-like equation. So a metric with exponentials along the diagonal, with a 1/R in the potential, might work.  Here is the exponential metric:

metric=(({{(E^(((((Sqrt[G]) q - 
                            G M))/((c^2) (Sqrt[(x^2) + (y^2) + (z^2)]))))), 0, 0, 0}, {0, -(E^(((2 (-(Sqrt[G]) q - G M))/((c^2) (Sqrt[(x^2) + (y^2) + (z^2)]))))), 0, 0}, {0, 0, -(E^(((2 (-(Sqrt[G]) q - G M))/((c^2) (Sqrt[(x^2) + (y^2) + (z^2)]))))), 0}, {0, 0, 0, -(E^(((2 (-(Sqrt[G]) q - G M))/((c^2) (Sqrt[(x^2) + (y^2) + (z^2)])))))}}));

g00 = metric[[1, 1]]

This is a singular solution:

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

If instead, the gauge we had chosen a flat, Euclidean metric, then the equation for the mass charge density would have been a Poisson equation, (Subscript[ρ, m]) = (∇^2) ϕ, which has the singular charge/R solution.

As an added check, use a tensor notebook written by Mathhew Headrick of MIT's Center for Theoretical Physics to confirm the calculation:

coord = {t, x, y, z};
<< diffgeo5.m
Sum[
      Simplify[D[Christoffel[[i, 1]][[1]], t] + 
          D[Christoffel[[i, 1]][[2]], x] + 
          D[Christoffel[[i, 1]][[3]], y] + 
          D[Christoffel[[i, 1]][[4]], z]], {i, 4}]

I like to confirm the result was not trivial (meaning all the components were not zero):

Table[Christoffel[[i, 1]][[1]], {i, 4}]

[A personal note: I had first derived the exponential metric in the way outlined in the next part of this notebook. Then a mere four years later, realized that the divergence of the Christoffel symbol had to return the mass density. Having never calculated a Christoffel symbol for any metric, all of my work for four years was on the line, in an equation I did not know how to calculate! If it did not work out, then I would have to tell everyone this line of research was wrong. Needless to say, I was relieved when the result turned out correctly.]

The best known show-stopping problem with a rank one field equation originates from solutions to the 4D wave equation, (Jq^μ) - (Jm^μ) = (□^2) (A^μ). A class of solutions has the form of an inverse distance squared.  To generate a force from such a potential, take the derivative.  The resulting force law has an inverse distance cubed dependence, so this potential is obviously not physical. Newtonian gravity is not an inverse cube force law, never was, never will be.

It has been shown what conditions generate Newton's field equations, (Jm^0) = (∇^2) (A^0). A potential that solves this equation has an inverse distance dependence, so the spatial derivative will be an inverse square, the correct form needed for a gravitational force law.  Imagine that (∂^2) (ϕ/∂)(t^2) is an incredibly small, but non-zero number.  There should be a smooth transition from the zero to non-zero situation, not a dramatic breakdown from inverse distance squared to an inverse cubed force law. Gravity is far weaker than the other three known forces.  Spacetime is hardly curved at all by the mass charge around us.  This suggests perturbation theory should be applied to the problem.  Given the way the Earth has wobbled around the Sun for four billion years, that kind of motion suggests a simple harmonic oscillator of some sort.  The potential must solve the field equations.  The derivative of the potential must under classical conditions have a (M/(R^2)) dependence. One possibility would involve a simple harmonic oscillator, with a spring constant related directly to the source mass (k = ((G M)/(c^2)) has units of L), over a distance squared (units of L^2).  If so, the derivative of the potential would have units of inverse distance.  A dimensionless potential would have a spatial derivative with units of inverse distance.  The hunt is on for a dimensionless perturbation that solves the field equations.  A historical note: when people first worked with the 4D wave equation in the eighteen hundreds, they would not have considered the notion of geometric length for a mass, an idea that arose from general relativity.

Let's try and keep this as simple as possible, without being too simple.  The idea is to study something a small step away from classical Newtonian gravitational physics - neutral, spherically symmetric, not rotating, a 1/R potential - with a modification to include a small contribution from time (translation: |R| >>> |ct|, but ct is not zero).  One such potential that solves the vacuum field equations is:

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

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[a1]

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

This is a well known result for 4D waves - an inverse distance squared solves the equation instead of an inverse distance like the 3D Poisson equation.  It is interesting that Planck's constant h is required to get the units correct since the presence of this constant is a sign of the domain of quantum mechanics.  The goal is to study a potential with the form of a normalized simple harmonic oscillator.  Do these steps together: create a small linear perturbation, and normalize the potential to the non-perturbed distance. [(σ^2) is the Lorentz-invariant distance interval, (x^2) + (y^2) + (z^2) - (c^2) (t^2)]:

a2 = {(((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[a2]

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

The potential a2 is a perturbation normalized to itself, using different constants to keep the units consistent.  Planck's constant is no longer required, indicating the problem has gone back to the a relativistic gravity domain (there is c and G).  The potential contains a shift (+1) and a rescaling ((k/(σ^2))), neither of which affect the solution to the differential equation.  Take the derivative of a2 with respect to time and space to first order in the spring constant k:

Expand[D[a2, 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[a2, t]] /. sublist

If the spring constant was equal to the geometric length of the source mass (k = ((G M)/(c^2)) which has units of distance), then:

Expand[D[a2, t]] /. sublist /. 
      k -> ((G M)/(c^2))
Expand[D[a2, x]] /. sublist /. 
      k -> ((G M)/(c^2))
Expand[D[a2, y]] /. sublist /. 
      k -> ((G M)/(c^2))
Expand[D[a2, z]] /. sublist /. k -> ((G M)/(c^2))

In the classical limit, ct -> 0 so σ -> R.  This has the form sought: the derivative of the potential is the source mass over a distance squared, so this could be a physically relevant potential.

Now we might understand why brilliant people in the nineteenth century passed by this possibility.  Perturbation theory is not considered "elegant".  It is merely practical.  To associate a source mass M with a spring constant with units of distance (((G M)/(c^2))) was not in the lexicon of the conceivable until the geometric approach of general relativity.

Even if researchers today were to go through that much analysis to get the right inverse distance dependence for the derivative of the potential, they would not be satisfied because experimentalists have confirmed small changes are required for Newton's field equation due to the finite speed of gravity.  What I hope to find is a potential, when put in a force equation, can generate a metric equation consistent with all weak gravitational field tests of general relativity.  From the modern viewpoint, this goal is odd, because general relativity tells us there is no force of gravity, only a dynamic change in the metric of spacetime.  I hope to demonstrate a force equation whose solution for a particular potential is a change in the metric of spacetime. That is what the algebra will do.

The first problem is to create a potential that only involves gravity, no EM.  A way to do this is with a potential whose derivatives lie along the diagonal of the field strength tensor:

(Subscript[∇, μ]) (A^ν)=(({{f, 0, 0, 0}, {0, g, 0, 0}, {0, 0, h, 0}, {0, 0, 0, j}}))
where f, g, h, 
    j are the appropriate derivatives of (A^ν).

The normalized perturbation potential worked with earlier would have off-diagonal elements, and thus represent a system with mass and electric fields.  The task is to work only with gravity.  The algebraic question becomes how to create a potential whose derivatives only lie on the diagonal, nothing off-diagonal.  Only first order terms of the spring constant k are being kept for the derivative. Changing the sign of the spring constant k does not affect the solution of the 4D wave differential equations, but does change the derivative of the potential.  For example, the following sum of two normalized perturbation potentials solves the 4D wave equation:

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]

MatrixForm[
      Expand[Simplify[contraD[diagonalSHO] /. sublist /. 
            k -> ((G M)/(c^2))]]]

Only the terms to first order in M will make a significant contribution:

MatrixForm[
      Normal[Series[
          Expand[Simplify[contraD[diagonalSHO] /. sublist /. 
                k -> ((G M)/(c^2))]], {M, 0, 1}]]]

It is ironic how much work is required to generate the identity matrix times a common factor.

The force equation for mass and electric charge will have the same form as the Lorentz force equation for electric charge alone:

Lorentz force equation for EM
\(F = \(\(1\/c\) 
          q\ \(U\_ν\) \((∂\^μ A\^ν - ∂\^ν A\^μ)\) = \(∂\ m\ U\^μ\)\/∂\ τ\ \ \ \ \ \ \ \ where\ U\^ν\ is\ the\ relativistic\ 4 - 
          velocity\)\)
Lorentz force equation for GEM
\(F = \(\(1\/c\) \((q - \(\@G\) 
                m)\)\ \(U\_ν\) ∇\^μ A\^ν = \(∂\ m\ U\^μ\)\/∂\ τ\)\)

Assume that there is no electric charge (q=0).  Plug in the derivative of the normalized potential which applies only to gravity from above:

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[v, 
          Expand[contraD[diagonalSHO]] /. sublist]]

Note G has vanished.  Don't worry, it is in the spring constant k. The relativistic force F is written in terms of a derivative with respect to the interval τ.  τ is related to σ by the equation (((cτ))^2) = -(σ^2).  Spherical symmetry will be presumed from here on, so (U_1, U_2, U_3) will be written as U_R.  Make these changes:

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

Apply the chain rule to the change in momentum.

D[m[τ] (U^ν)[τ], τ]

D[m[τ] (Subscript[U, μ])[τ], τ]

In this section, we assume that test mass does not change in spacetime. To find a new constant velocity solution in a later section, this assumption will not be made.

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

Presume the equivalence principle so the gravitational test mass m on the left side of the force equation cancels the inertial mass m on the right.  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])[τ]}, τ]

The next task is to eliminate the constants of integration. Contract the relativistic velocity 4-vector solution.

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

In flat spacetime, k->0

Limit[(E^(-((2 k)/(c τ)))) C[1]^2 - (E^(((2 k)/(c τ)))) C[2]^2, 
      k -> 0]

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

flatSpacetimeConstraint = { C[1] -> c (dt/dtau), 
          C[2] -> (dR/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

Since perturbation theory has been used, it has already been assumed that the metric is close to flat, thus the contracted velocity is approximately (c^2). 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))

The magnitude of tau will be the same as sigma (|c tau| = |sigma|~=R). To make the metric real, choose sigma to be imaginary, so that cτ ≂ R.  Make the substitutions τ -> (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)

The equation has the same form as the exponential metric discussed earlier

MatrixForm[
      DiagonalMatrix[{(E^(-((2 G M)/((c^2) R)))), - (E^(((2 G M)/((c^2) R)))), - (E^(((2 G M)/((c^2) R)))), - (E^(((2 G M)/((c^2) R))))}]]

\(Take\ the\ Taylor\ series\ expansion\ for\ a\ small\ value\ of\ the\ gravitational\ field\ \(\(\(G\ M\)\/\(c\^2\ R\)\)\(.\)\)\)

(Series[
            (dt^2) (E^(-2 ϕ)) - 
              ((dR^2)/(c^2)) (E^(2 ϕ)), {ϕ, 0, 
              2}] /. ϕ -> ((G M)/((c^2) R))) == (dtau^2)

Rewrite.

(1 - 2 (( G M)/((c^2) R)) + 2 (((((G M)/((c^2) R))))^2)) 
            (dt^2) - (1 + 2 (( G M)/((c^2) R))) ((dR^2)/(c^2)) == 
        (dtau^2);

The Taylor series expansion of the isotropic Schwarzschild metric (MTW, exercise 31.7)

(Series[(((((1 - (ϕ/2))/(1 + (ϕ/2)))))^2) 
                (dt^2) - (((1 + (ϕ/2)))^4) ((dR^2)/(c^2)), {ϕ, 0, 
              2}] /. ϕ -> ((G M)/((c^2) R))) == (dtau^2)

Rewrite.

(1 - 2 (( G M)/((c^2) R)) + 2 (((((G M)/((c^2) R))))^2)) 
            (dt^2) - (1 + 2 (( G M)/((c^2) R))) ((dR^2)/(c^2)) == 
        (dtau^2);

This is exactly the terms tested and confirmed for weak field tests of general relativity, so this proposal will pass all those tests. This is known as first-order Parametrized Post-Newtonian (1st order PPN) accuracy. I spoke to Clifford Will who said that tests to second order PPN accuracy are not being conducted or even planned. So what is the predicted difference at second order PPN accuracy? Reuben Epstein and Irwin Shapiro provide a formula (Phys. Rev. D, 22:2947, 1980):

\(Δθ = π \((2\  + 
            2  γ - β + \(3\/4\) ϵ)\) \(G\^2\ M\^2\)\/\(\(c\^4\) R\^2\)\)

For general relativity, γ, β and ϵ are all equal to one:

Δθ = π (2  + 
              2  γ - β + (3/4) ϵ) (((G^2) (M^2))/((c^4) (R^2))) 
          microarcseconds /. {γ -> 1, β -> 
            1, ϵ -> 1, M -> 1.98 (10^30) kg, 
          R -> 6.96 (10^8) m, 
          G -> 6.67 (10^(-11)) ((m^3)/(kg (s^2))), 
          c -> 3 (10^8) (m/s), 
          microarcseconds -> 
            1/(((2 π)/360) (1/60) (1/60) 
                  (10^(-6)))}

General relativity predicts bending at the solar limb at 10.8 microarcseconds. The bending for the exponential metric sets ϵ equal to 4/3:

Δθ = π (2  + 
              2  γ - β + (3/4) ϵ) (((G^2) (M^2))/((c^4) (R^2))) 
          microarcseconds /. {γ -> 1, β -> 
            1, ϵ -> 4/3, M -> 1.98 (10^30) kg, 
          R -> 6.96 (10^8) m, 
          G -> 6.67 (10^(-11)) ((m^3)/(kg (s^2))), 
          c -> 3 (10^8) (m/s), 
          microarcseconds -> 
            1/(((2 π)/360) (1/60) (1/60) 
                  (10^(-6)))}

The exponential unified field theory metric predicts 0.7 microarcseconds more bending than the Schwarzschild metric of general relativity.The current measurements of light bending are on the order of a hundred microarcseconds. At the microsecond level, other effects like the speed of the rotation of the Sun and its quadrupole moment come into play. Still this is an important milestone: the exponential metric is not a solution to the Einstein field equations, but could be confirmed or rejected on experimental grounds.

V: Classical constant velocity solution

We first need to derive Newton's gravitational force law from this completely relativistic one.

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

Newton's classical force law is conservative, so the first terms of the above equation are zero. Presume the change in mass term contributes nothing. The spring constant k becomes the gravitational length of the source mass. The distance sigma becomes R.

conservativeForce = {(Subscript[U, 0])[τ] -> 
            0, (Subscript[U, 0])'[τ] -> 0};
noMassChange = {m'[τ] -> 0,  
          m[τ] -> m};
springIsSourceMass = {k -> ((G M)/(c^2))};
sigmaToR = {σ -> R};

{-((c k m[τ] (Subscript[U, 0])[τ])/(σ^2)), ((c k m[τ] (Subscript[U, R])[τ])/(σ^2))} == {D[
                  m[τ] (Subscript[U, 0])[τ], τ], 
                D[m[τ] (Subscript[U, R])[τ], τ]} /. noMassChange /. 
          conservativeForce /. springIsSourceMass /. sigmaToR

We now need to break spacetime symmetry. We can no longer use a relativistic 4-velocity or 4-acceleration. The question is what is now the appropriate derivatives and directions for those derivatives? Newton's law describes a static force field, so the interval tau has the same magnitude as the absolute value of the distance, |R|.

{D[t[τ], τ], 
        c D[R[τ], τ]} /. {τ -> R}

In classical physics, time is independent of space, so the gamma term here, ∂(t/∂)R, is zero. The other term is a unit vector in the R direction. This says that change only happens along the direction of R, a reasonable statement.

{D[t[τ], τ], 
          (c^2) D[R[τ], {τ, 2}]} /. {τ -> R} /. t'[
          R] -> 0

This acceleration still is not classical because it contains the constant (c^2). One way to eliminate it is to substitute R/c->it. Do that twice, and in pops a minus sign, out go the c's.

{D[t[τ], τ], 
            (c^2) D[R[τ], {τ, 2}]} /. {τ -> 
              R} /. t'[R] -> 0 /. R''[
          R] -> -(1/(c^2)) R''[t]

This is the classical acceleration. Plug this substitutions into the relativistic force law:

\(\(changeOnlyAlongRhat = {U\_R[τ] -> 
            c\ R\&~};\)\)
velocity2dRdtau = {(Subscript[U, R]) -> R'[τ]};
dtau2t = {R''[τ] -> -R''[t]};

{-((c k m[τ] (Subscript[U, 0])[τ])/(σ^2)), ((c k m[τ] (Subscript[U, R])[τ])/(σ^2))} == {D[
                        m[τ] (Subscript[U, 0])[τ], τ], 
                      D[m[τ] (Subscript[U, R])[τ], τ]} /. noMassChange /. 
                conservativeForce /. springIsSourceMass /. sigmaToR /. 
          changeOnlyAlongRhat /. velocity2dRdtau /. dtau2t

This is Newton's gravitational force law.

One interesting thing to explore is the consequences of not assuming the distribution of mass in fact does change in spacetime. The derivation relied on assume a static system, so the change in mass with respect to time must be zero. Since the derivative is with respect to spacetime, the effect of gravity could also be on the distribution of mass in space. Excluding point sources, it is clear that the density of mass does vary with respect to distance.

If there is a new effect of the gravitational force with respect to distance, the vector must point along the velocity vector, not Rhat.

\(\({0, \(G\ m\ M  R\&~\)\/R\^2} == {0, \(-D[m[t] D[R[t], t], 
                t]\)} /. R\&~ -> R\&~ + V\&~\) /. \(m'\)[t] -> 
        c\ \(m'\)[R]\)

DSolve[((G M m[R])/(c V (R^2))) == -D[m[R], R], m[R], 
      R]

The c in the new term is important because it says the new effect will be small (most of any effect would be accounted for by c). This small effect will have a 1/R dependence, not 1/(R^2). This is precisely the kind of form of an equation that the Modification of Newtonian Dynamics (MOND): a small factor that eventually makes a big contribution because of the 1/R dependence.

So far, I have been unable to apply this new solution to a problem in astronomy.  The math is a bit tricky, and my skills are limited. It is absurdly rare to find a new class of solutions involving gravity, but it must be applied to a known situation before the approach has any possibility of being accepted.

VI: The Stress-Energy Tensor

The stress-energy tensor is a second rank symmetric tensor. Here is formula to calculate it:

\(T\^μν = \(∂ℒ\/∂\(∇\_μ A\_σ\)\) ∇\^ν A\_σ - \(\(g\^μν\) \(ℒ\)\(\ \)\)\)

Calculate the derivative of the Lagrangian:

\(\(∂ℒ\/∂\(∇\_μ A\_σ\)\) ∇\^ν A\_σ = ∇\^μ A\^σ ∇\^ν A\_σ\)

The two tensors that contract are easy enough to calculate separately:

MatrixForm[contraD[A]]

MatrixForm[contraD[covariant[A]]]

Make a function to do the right kind of contraction.

guv=(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}));

MinkowskiMetric=(({{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}}));

contractMM[contraD[A], coD[A]]

contractMM2M[m1_, m2_] := 
      Module[{u}, 
        Table[{m1[[u, 1]] m2[[1, 1]] - 
              m1[[u, 2]] m2[[1, 2]] - m1[[u, 3]] m2[[1, 3]] - 
              m1[[u, 4]] m2[[1, 4]], -m1[[u, 1]] 
                m2[[2, 1]] + m1[[u, 2]] m2[[2, 2]] + 
              m1[[u, 3]] m2[[2, 3]] + 
              m1[[u, 4]] m2[[2, 4]], -m1[[u, 1]] 
                m2[[3, 1]] + m1[[u, 2]] m2[[3, 2]] + 
              m1[[u, 3]] m2[[3, 3]] + 
              m1[[u, 4]] m2[[3, 4]], -m1[[u, 1]] 
                m2[[4, 1]] + m1[[u, 2]] m2[[4, 2]] + 
              m1[[u, 3]] m2[[4, 3]] + 
              m1[[u, 4]] m2[[4, 4]]}, {u, 4}]]

MatrixForm[
      Simplify[Expand[-(1/(2  (c^2))) 
              contractMM2M[contraD[A], 
                contraD[covariantvec[A]]]  + (1/(2  (c^2))) 
              guv contractMM[contraD[A], coD[A]]]]]

Expand[%[[1, 1]]]

This is a symmetric matrix, as it must be to conserve angular momentum. In classical electrodynamics, the canonical definition of the stress tensor density does not generate a symmetric matrix, and this must be corrected.

Write the first term using the fields:

Expand[(1/(c^2)) ((1/2) g[[1]]^2 - (1/2) 
              g . g - (1/8) (Efield + e) . (Efield + e) - (1/4) 
              B . B  - (1/4) b . b)]

For classical EM, the expected result is the sum of the squares of the electric and magnetic fields.  I am not sure what to make of all the differences in this expression.

VII: Quantization

In this section, I will see if it is reasonable to hope that the modes of radiation for the field equations can be quantized. What is one way to tell if modes of radiation for a set of field equations generated from a Lagrange density cannot be quantized? Classically, energy and momentum are numbers.  In quantum mechanics, energy and momentum are no longer numbers.  Instead they are operators which act on a wave function.  The average value of energy and momentum can be measured.  An operator cannot commute with its complementary operator.  For example, the position operator x and the momentum operator Px, their commutator, [x, Px] = (x Px - Px x), cannot equal zero.  If the operator Px was always zero, the commutator [x, Px] would also be zero.  For such a situation, the proposal could not be quantized.

Let's examine the situation for the classical EM Lagrange density:

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

Define a function to calculate the generalized momentum density:  The only tricky part here is that Mathematica does not want to treat a partial derivative as a derivative variable, so substitutions need to be made, then undone.

momentum[L_] := 
      Module[{noPd, toPd, mo}, 
        noPd = {(∂_t) ϕ[t, x, y, z] -> 
              dphidt, (∂_x) ϕ[t, x, y, z] -> 
              dphidx, (∂_y) ϕ[t, x, y, z] -> 
              dphidy, (∂_z) ϕ[t, x, y, z] -> 
              dphidz, (∂_t) Ax[t, x, y, z] -> 
              dAxdt, (∂_x) Ax[t, x, y, z] -> 
              dAxdx, (∂_y) Ax[t, x, y, z] -> 
              dAxdy, (∂_z) Ax[t, x, y, z] -> 
              dAxdz, (∂_t) Ay[t, x, y, z] -> 
              dAydt, (∂_x) Ay[t, x, y, z] -> 
              dAydx, (∂_y) Ay[t, x, y, z] -> 
              dAydy, (∂_z) Ay[t, x, y, z] -> 
              dAydz, (∂_t) Az[t, x, y, z] -> 
              dAzdt, (∂_x) Az[t, x, y, z] -> 
              dAzdx, (∂_y) Az[t, x, y, z] -> 
              dAzdy, (∂_z) Az[t, x, y, z] -> 
              dAzdz}; 
toPd = {dphidt -> (∂_t) ϕ[t, x, y, z], 
            dphidx -> (∂_x) ϕ[t, x, y, z], 
            dphidy -> (∂_y) ϕ[t, x, y, z], 
            dphidz -> (∂_z) ϕ[t, x, y, z], 
            dAxdt -> (∂_t) Ax[t, x, y, z], 
            dAxdx -> (∂_x) Ax[t, x, y, z], 
            dAxdy -> (∂_y) Ax[t, x, y, z], 
            dAxdz -> (∂_z) Ax[t, x, y, z], 
            dAydt -> (∂_t) Ay[t, x, y, z], 
            dAydx -> (∂_x) Ay[t, x, y, z], 
            dAydy -> (∂_y) Ay[t, x, y, z], 
            dAydz -> (∂_z) Ay[t, x, y, z], 
            dAzdt -> (∂_t) Az[t, x, y, z], 
            dAzdx -> (∂_x) Az[t, x, y, z], 
            dAzdy -> (∂_y) Az[t, x, y, z], 
            dAzdz -> (∂_z) Az[t, x, y, 
                  z]}; 
mo = {D[L /. noPd, dphidt], 
            D[L /. noPd, dAxdt], D[L /. noPd, dAydt], 
            D[L /. noPd, dAzdt]}; 
mo /. 
          toPd
]

Calculate the generalized 4-momentum of the classical EM Lagrange density.

momentum[LEM]

The modes of emission cannot be quantized due to that zero. One way around this problem is to fix the gauge.  One covariant approach is to use the Gupta/Bleuler Lagrange density:

Expand[(LGB = -covariant[(Subscript[J, q])] . A/
                c - (1/(2 
                       (c^2))) (((D[A[[1]], t] + c D[A[[2]], x] + c D[A[[3]], y] + c D[A[[4]], z]))^2) - 
              Expand[(1/(4 (c^2))) (contractMM[
                      contraD[A] - contraDvu[A], 
                      coD[covariant[A]] - 
                        coDvu[covariant[A]]])]) /. {ϕ[t, x, y, 
              z] -> ϕ, Ax[t, x, y, z] -> Ax, 
          Ay[t, x, y, z] -> Ay, Az[t, x, y, z] -> Az}]

Calculate its generalized 4-momentum.

momentum[LGB]

These field equations can be quantized, but run into a different technical problem.  The field strength tensor is second rank and antisymmetric, so will be represented by a spin one field where like charges repel. The scalar mode of emission for a spin one field could have a negative energy density, and that makes no sense. An additional constraint is required to make the scalar and longitudinal spin 1 fields virtual.

Calculate the generalized 4-momentum of the GEM Lagrange density.

momentum[LGEM]

The GEM field equations do not have any zeros in the 4-momentum density, so it is possible to quantize the modes of emission.

The GEM field is NOT all a spin 1 field.  The field strength tensor has two parts.  The antisymmetric second rank field strength tensor will be represented by a spin 1 field for EM where like charges repel. The symmetric second rank field strength tensor will be represented by a spin 2 field for gravity where like charges attract.  These modes of emission are scalar and longitudinal. The scalar mode of a spin 2 field will not have the negative energy density problem. We know that gravity has a classical longitudinal wave behavior. If one drills a hole through the center of the Earth, and creates a vacuum in that tunnel, a ball dropped in the tunnel will oscillate with a period of about an hour and a half. The direction the particle accelerates is the direction it is moving, so the wave is longitudinal. Should we ever detect the polarization of gravity waves, this unified field proposal predicts the modes will either be scalar or longitudinal, and not transverse as is predicted by general relativity. As such, the polarization of gravity waves is a clear way to confirm or reject this proposal.