I will try to keep my Newtonian wig on. Digressions will be in well-defined sidebars. Time to serve the meat and potatoes.

The Action for the Newtonian Gravitational Field

Two things go into this action: a mass density and a scalar field:

\begin{align*} S_{_{NGF}}&= \int \mathcal{L}(\Phi, \nabla \Phi)\\ &= \int d^4 x \left(-\rho\, \Phi - \frac{c^2}{2 G} (\nabla \Phi)^2\right ) \end{align*}

I have chosen to work with a dimensionless source potential. Vary the action with respect to the potential. Use the Euler-Lagrange equation to generate the field equations:

\begin{align*} \frac{d \mathcal{L}}{d \Phi} &= \nabla \left( \frac{\partial \mathcal{L}}{\partial \nabla \Phi}\right) \\ \rho &= \frac{c^2}{G} \nabla^2 \Phi \end{align*}

The general solution is - wait, there is no general solution.  Solutions always depend on the boundary conditions.  Let's work on a simple case: outside of a spherically symmetric, non-rotating, uncharged, not-too-dense source. The potential in that situation is:

$\Phi = k_1 \frac{M}{R} + k_2$

[My bad here: the source mass M is part of the integration constant that gets set by the boundary conditions. I was too used to seeing M/R that I just put it in out of habit. LFB pointed out this error in the comments below.]

The integration constants cannot be pinned down with only one equation and two unknowns.

The Equations of Motion of a Tiny Test Mass

The tiny test mass is so small it does not materially change the field equations above.  The action is a function of a parameterized position and velocity:

\begin{align*} S_{_{NGEq~of~Motion}}&= \int \mathcal{L}\left(\lambda, R(\lambda), \frac{d R(\lambda)}{d t} \right ) \\ &= \int d\lambda \left(\frac{m}{2 c^2} \left( \frac{d R(\lambda)}{d t} \right )^2 - m \Phi(R(\lambda))\right ) \end{align*}

The speed of light is finite in Newtonian mechanics.  Both terms need to have the same units, so that is why the c2 was used.

The tiny test mass shows up in both terms: as the inertial mass and as the passive gravitational mass that couples to the gravity field.  The tiny test mass is acting like a charge in EM.  Vary the action with respect to position.  Use the Euler-Lagrange equation to generate the equations of motion:

\begin{align*} \frac{d}{d \lambda} \left( \frac{d \mathcal{L}}{d \left( \frac{d R}{d \lambda}\right)}\right) &= \frac{d \mathcal{L}}{d R} \\ \frac{m}{c^2} \frac{d^2 R}{d \lambda^2} &= -m \nabla \Phi \\ &= -m \nabla\left( k_1 \frac{1}{R} + k_2\right) \end{align*}

[Note: modified from the original which had k1 M/R, see note above.]

Gravity is an attractive force, so the acceleration should be negative. Along with the requirement that the potential is dimensionless, that sets k1=-GM/c2.  It is unclear to me how to fix the value of the second integration constant. It is a good thing to have a potential that is always positive. If the second constant was say 17, that would be the case. In a metric theory, the result is constrained to match the Minkowski metric without any source mass. The Newtonian wig prevents me from using that justification :-) To make the potential simple and positive, I will set k2 equal to one:

\begin{align*} \frac{m}{c^2}\frac{d^2 R}{d \lambda^2} &= -m \nabla \left(- \frac{G M}{c^2 R} + 1 \right )\\ m\frac{d^2 R}{d t^2} &=-m \frac{G M}{R^2} \end{align*}

The equation of motion for a tiny test particle around a spherically symmetric, uncharged, non-rotating, not-too-dense source look like the force equation Newton himself wrote so many years ago. To be honest, I have always thought of Newton's gravity force law as something chiseled in stone. Yet writing out all the requirements for the equations of motion makes the law feel less universal. This feels like a shift, from the one gravity force LAW, to a situation where one has to find out the equations to be used based on the conditions at hand and then calculate the motion of the system.

Sidebar: Generalizing the Source

In the action for the field equations, I could have used one term of the stress-energy tensor:

$\rho = T_{00}$

A lesson from general relativity is that all energy should be accounted for as a driver of gravity. It is sensible in terms of the math to use the square root of the contraction of the stress-energy tensor?

$\rho \rightarrow \sqrt{T_{\mu \nu}T^{\mu \nu}}$

The value of rho depends on a metric via the contraction of a 4-momentum. The value of the tensor contraction has two contractions that one then takes the square root. The difference between the two is the inclusion of the other terms in the stress-energy tensor. Is anyone familiar with experimental data for those contributions?

One thing I am unsure about is whether it will change how the new source transforms under a Lorentz transformation. The mass density rho will be a Lorentz invariant mass, and an inverse volume, so will change based on the volume element.  It looks to my eye that the contraction of the stress-energy tensor may do the same thing. Comments?

Note to commentors: please stick to the topic at hand, Newton's model for gravity written as an action. I reserve the right to delete comments and will do so after leaving my own explanation as to why.  If you have a great solution to a problem, congrats, blog about it on this website, just not in this particular blog.