Well, you see there's a problem there. The actual kinematic symmetry group for non-relativistic physics is not the Galilei group, but the Bargmann group - its central extension. This is best seen by how mass $m$, momentum $𝐩 = \left(p_x, p_y, p_z\right)$ and kinetic energy $H$ transform:
$$Δm = 0, \hspace 1em Δ𝐩 = 𝞈×𝐩 - 𝞄m, \hspace 1em ΔH = -𝞄·𝐩,$$
where $𝞈$ denotes an infinitesimal rotation and $𝞄$ an infinitesimal (Galilei) boost.
The problem is actually much bigger than I'll be letting on here and I'll only mention this briefly in passing: there is no real $c = ∞$ limit, because there's a total mismatch between the Galilei (and Bargmann) group versus the Poincaré group. The former has a "wild" representation type, which effectively means that its representations are far too many to even be classified. No matter what kind of $c = ∞$ limit you impose, there will always be some representations (meaning: almost all of them) that lie outside the scope of the reverse translation (i.e. the "relativization" process) and are lost. They get no "relativistic" version. In that sense, Relativity is a permanently incomplete paradigm. It can't capture and "relativize" the full complexity of what is present in the non-relativistic world - and some of what's lost may be physically relevant.
In here and below, it's easiest to discuss the transforms in infinitesimal form, since the algebra is linear and simpler than with transforms in finite form.
With Galilei transforms, you're forced to drop down to $m = 0$ and the transforms reduce to that for the 4-vector $(𝐩,H)$:
$$Δ𝐩 = 𝞈×𝐩, \hspace 1em ΔH = -𝞄·𝐩.$$
By comparison, the action on the coordinates, their differentials and the corresponding differential operators is:
$$
Δ𝐫 = 𝞈×𝐫 - 𝞄t + 𝝴, \hspace 1em Δt = τ, \\
Δ(d𝐫) = 𝞈×(d𝐫) - 𝞄(dt), \hspace 1em Δ(dt) = 0, \\
Δ(∇) = 𝞈×∇, \hspace 1em Δ\left(\frac{∂}{∂t}\right) = 𝞄·∇
$$
the last of these being derived from the second of these by the requirement that
$$d𝐫·∇ + dt \frac{∂}{∂t}$$
be invariant. For the coordinates, you also have the infinitesimal spatial translations $𝝴$ and infinitesimal time translations $τ$. That yields the inhomogeneous extension of the Galilei group (which is normally what is known as the Galilei group).
Thus, the 4-vector $(𝐩,-H)$ - when $m = 0$ - transforms as does $(∇, ∂/∂t)$. In turn, this means $𝐩·d𝐫 - H dt$ is a Galilei invariant.
But, when $m ≠ 0$, you have a mismatch: there are 5 components in $(m,𝐩,H)$, but only 4 in $(d𝐫,dt)$ and $(∇, ∂/∂t)$. What does $m$ go with? And - more importantly - how do you recover the action on $m$ in coordinate form?
The best way to answer that is to turn $m$ back on, i.e. let $m ≠ 0$, and see what happens with the (now, former) invariant $𝐩·d𝐫 - H dt$:
$$\begin{align}
Δ(𝐩·d𝐫 - H dt) &= (Δ𝐩)·d𝐫 + 𝐩·Δ(d𝐫) - (ΔH) dt - H (dt) \\
&= (𝞈×𝐩 - 𝞄m)·d𝐫 + 𝐩·(𝞈×(d𝐫) - 𝞄(dt)) - (-𝞄·𝐩) dt - H (0) \\
&= 𝞈×𝐩·d𝐫 + 𝐩·𝞈×d𝐫 - 𝞄m·d𝐫 - 𝐩·𝞄dt + 𝞄·𝐩dt.
\end{align}$$
After noting the cancellation
$$𝞈×𝐩·d𝐫 + 𝐩·𝞈×d𝐫 = 𝞈×𝐩·d𝐫 - 𝞈×𝐩·d𝐫 = 0,$$
this reduces to:
$$Δ(𝐩·d𝐫 - H dt) = m(-𝞄·d𝐫).$$
The answer of what $m$ goes with is something, $u$, that ought to transform as
$$Δu = 𝞄·𝐫 + ψ, \hspace 1em Δ(du) = 𝞄·d𝐫,$$
now with an additional transform $ψ$ for an infinitesimal $u$-translation.
Once you have that in place, you can create an invariant: $𝐩·d𝐫 - H dt + m du$. The inclusion of $u$ alters the transform for the differential operators, since we now require that this be the invariant:
$$d𝐫·∇ + dt \frac{∂}{∂t} + du \frac{∂}{∂u}.$$
Correspondingly, what we get is:
$$Δ(∇) = 𝞈×∇ - 𝞄\left(\frac{∂}{∂u}\right), \hspace 1em Δ\left(\frac{∂}{∂t}\right) = 𝞄·∇, \hspace 1em Δ\left(\frac{∂}{∂u}\right) = 0.$$
This is action of the Bargmann group, and the geometry that it goes with is not 4-dimensional at all, but 5-dimensional: the Bargmann geometry.
So, if you want a Newtonian analogue of General Relativity, cast in the same mold as General Relativity (and that's not just "Newton-Cartan geometry", which was first developed before people became aware that the Galilei group was the wrong symmetry group), then you have to make a few changes, so that an actual correspondence between General Relativity and its Newtonian version can be cleanly bridged.
The change actually has to be made on the General Relativity side, more so than on the Newtonian side; and I don't think that matter has been fully resolved. Just as Galilei was lifted to Bargmann, on the other side of the Correspondence Principle, Minkowski should also be lifted to the (unnamed) relativistic version of Bargmann; and you work from there.
Compare the situation with Relativity. First, note that the above transforms yield the following geometric invariants, in addition to the ones already cited above for the differential operator and the mass-momentum-energy 1-form:
$$
|d𝐫|^2 + 2 dt du, \hspace 1em dt, \\
∇^2 + 2\frac{∂}{∂t}\frac{∂}{∂u}, \hspace 1em \frac{∂}{∂u}, \\
p^2 - 2mH, \hspace 1em m.
$$
Is there anything like that in Relativity?
Well, let's look at the transform behavior of momentum and energy. Historically, instead of considering "total energy" $E$, one considered the "moving mass" $M$ and the kinetic energy $H$ the same as above, as well as the "rest mass" $m$. They have the following relations:
$$
M = αE, \hspace 1em H = E - m/α, \\
E = M/α, \hspace 1em m = M - αH.
$$
where I will use the parameter $α = (1/c)^2$ to distinguish between the relativistic and non-relativistic cases, the latter obtained (where possible) by setting $α = 0$.
The historical quantities $M$ and $H$ gave way to $m$ and $E$, over time - which obscured the connection between the relativistic and non-relativistic cases. What distinguished relativity is that it stipulated a non-trivial boost transform on $t$ and (going along with this) a non-trivial transform on $M$. For the coordinates and their differentials, the following was posed:
$$
Δ𝐫 = 𝞈×𝐫 - 𝞄t + 𝝴, \hspace 1em Δt = -α𝞄·𝐫 + τ, \\
Δ(d𝐫) = 𝞈×(d𝐫) - 𝞄(dt), \hspace 1em Δ(dt) = -α𝞄·(d𝐫), \\
Δ(∇) = 𝞈×∇ + α\left(\frac{∂}{∂t}\right), \hspace 1em Δ\left(\frac{∂}{∂t}\right) = 𝞄·∇,
$$
again, the action on the differential operators being determined by the requirement that $d𝐫·∇ + dt(∂/∂t)$ be invariant.
Unlike the Galilei case, we have an invariant with the momentum and energy - provided we take the "total energy" instead of the kinetic energy: $𝐩·d𝐫 - E dt$. The total energy and momentum, together as a 4-vector $(𝐩,-E)$, transforms as does $(∇, ∂/∂t)$, with the action:
$$Δ𝐩 = 𝞈×𝐩 - α𝞄E, \hspace 1em ΔE = -𝞄·𝐩.$$
Unlike the Galilei case, this includes a non-zero mass ... because of That Equation $E = Mc^2$, which we've been writing in disguised form as $M = αE$. So, we could just as easily write:
$$Δ𝐩 = 𝞈×𝐩 - 𝞄M, \hspace 1em ΔM = -α𝞄·𝐩,$$
and this shows that the "simultaneity-shift" $-α𝞄·𝐫$ over in $Δt$ found its way into the transform for the "moving mass" $M$. This is what lies behind mass-energy unification.
Because of that unification, what was 5: $(M=m,𝐩,H)$ became 4: $(𝐩,E)$.
We're going to revert and make it 5, again assuming $m$ is an invariant. The corresponding transforms then become:
$$Δ(M,𝐩,H) = (-α𝞄·𝐩, 𝞈×𝐩 - 𝞄M, -𝞄·𝐩),$$
and we have the following invariants:
$$|𝐩|^2 - 2MH + αH^2, \hspace 1em M - αH.$$
The former is 0, if we set the latter to $m$.
If we continue to impose the requirement that $𝐩·d𝐫 - H dt + m du$ remain invariant, then we get the (unnamed) relativistic version of Bargmann geometry, with the following action on the coordinates and their differential operators:
$$
Δ𝐫 = 𝞈×𝐫 - 𝞄t + 𝝴, \hspace 1em Δt = -α𝞄·𝐫 + τ, \hspace 1em Δu = 𝞄·𝐫, \\
Δ(d𝐫) = 𝞈×(d𝐫) - 𝞄(dt), \hspace 1em Δ(dt) = -α𝞄·(d𝐫), \hspace 1em Δ(du) = 𝞄·(d𝐫), \\
Δ(∇) = 𝞈×∇ + 𝞄\left(α\frac{∂}{∂t} - \frac{∂}{∂u}\right), \hspace 1em Δ\left(\frac{∂}{∂t}\right) = 𝞄·∇, \hspace 1em Δ\left(\frac{∂}{∂u}\right) = 0.
$$
The corresponding geometric invariants are:
$$
|d𝐫|^2 + 2 dt du + α du^2, \hspace 1em ds = dt + α du, \\
∇^2 + 2\frac{∂}{∂t}\frac{∂}{∂u} - α\left(\frac{∂}{∂t}\right)^2, \hspace 1em \frac{∂}{∂u}, \\
d𝐫·∇ + dt \frac{∂}{∂t} + du \frac{∂}{∂u}.
$$
Now, you can bridge the gap between the two: between the non-relativistic and relativistic cases.
What is the invariant $ds$ and what does $s = t + αu$ correspond to? In the non-relativistic case it $s = t$ is just the time. But, it's invariant. Therefore, it's a holdover of the "absolute time" of non-relativistic theory.
If you equate it with proper time and substitute into the Minkowski metric - written as a proper time metric:
$$-c^2 ds^2 = dx^2 + dy^2 + dz^2 - c^2 dt^2,$$
you get
$$dx^2 + dy^2 + dz^2 + 2 dt du + α du^2 = 0.$$
So, everything in 4D Minkowski space - all the dynamics in that space - are now taking place on light cones in the "relativistic" Bargmann geometry.
The difference between it and the Bargmann geometry can be seen by using the quadratic line element, itself, as a metric (the geometry is then 4+1 dimensional, regardless of what $α$ is), and apply it to the invariants $ds$ and $∂/∂u$. For $α = 0$, they are both null vectors in 5D. For $α > 0$, they are both space-like. In each case, the geometries are picking out an invariant direction and are threaded by an invariant vector field. Had the case instead been $α < 0$, it would have been a time-like invariant direction and vector field, and the reduced 4D geometry would have just 4D Euclidean space.
As an example, the main example that all of this runs up to, consider the Schwarzschild solution, now written as a proper time line element with $s = t + αu$ being the proper time:
$$
-\frac{ds^2}{α} = -\left(1 + 2V\right)\frac{dt^2}{α} + \frac{dr^2}{1 + 2αV} + r^2 \left((dθ)^2 + (\sin{θ} dφ)^2\right),
$$
which we'll write in a more Cartesian form by replacing the angular coordinates:
$$r^2 \left((dθ)^2 + (\sin{θ} dφ)^2\right) = dx^2 + dy^2 + dz^2 - dr^2,$$
as
$$
-\frac{ds^2}{α} = -\left(1 + 2αV\right)\frac{dt^2}{α} + dr^2\left(\frac{1}{1 + 2αV} - 1\right) + dx^2 + dy^2 + dz^2.
$$
Here,
$$V = -\frac{GM}{r}$$
is the gravitational potential per unit mass for a body centered on $r = 0$, with mass $M$.
After substituting for $ds$ and making the following simplification
$$\frac{1}{1 + 2αV} - 1 = -\frac{2αV}{1 + 2αV},$$
we get the following line element / constraint:
$$dx^2 + dy^2 + dz^2 + 2 dt du + α du^2 - 2V dt^2 - \frac{2αV dr^2}{1 + 2αV} = 0.$$
The non-relativistic version of this, obtained by setting $α = 0$, is ... The Line Element For Newtonian Gravity:
$$dx^2 + dy^2 + dz^2 + 2 dt du - 2V dt^2 = 0.$$
In the general case, you're essentially talking about an "Einstein-Aether" Theory in a 4+1 dimensional ambient geometry: a theory posed on a geometry that has an invariant field. For the non-relativistic case, this field is a null field.
You can, indeed, write out the Einstein-Hilbert action in 4+1 dimensions and apply the condition that the $ds$ direction and $∂/∂u$ be invariant. It is enough to give you Newtonian gravity, but you lose the connection that Einstein's Equations have to fluid dynamics in the $α = 0$ limit.
The 5D method has been done up to a certain point.
Bargmann structures and Newton-Cartan theory
C. Duval, G. Burdet, H. P. Künzle, and M. Perrin
Phys. Rev. D 31, 1841 – Published 15 April 1985
https://journals.aps.org/prd/abstract/10.1103/PhysRevD.31.1841
but there, too, the authors point out the obstruction.
"VI. Newton's Field Equations"
"Strangely enough, Newton's field equations" [their (1.15) and (1.17) cited] "cannot be easily derived from a specific space-time Lagrangian density. It seems that there might exist some puzzling geometric obstruction to the existence of a well-defined variatonal problem in the four-dimensional picture. Also the fact that only matter density enters the source term of Newton's field equations" [i.e. that the integration of fluid dynamics with the field equations, that Einstein's equations has, is lost] "has not been quite understood so far in a covariant formalism. The role of the mass flow and stress-energy tensor of matter distributions still remains to be clarified."
This is really an allusion to what I pointed out up front: stuff gets lost in the correspondence limit - and its inverse. When we introduce $u$, we recovered a relativistic version of the Newtonian absolute time. This also recovered a Newtonian version of "time dilation". The $u$ coordinate is just the time dilation $s - t$, itself, scaled up by $c^2$ so that it survives the $c = ∞$ limit.
It seems that there are more residual terms at the $α^2$ order that need to be taken into consideration, as well. The $du$ coordinate is the $α^1$ residue of $dt$, and $H$ the $α^1$ residue of $M$. One way might be to add in similar residues for the momentum and spatial coordinates.
An attempt to gauge the Galilei/Bargmann group, bringing in more residuals, may be found here:
Newtonian Gravitation as a Gauge Theory.
https://arxiv.org/abs/gr-qc/9405046
Generalized Newtonian Theories Of Gravitation
https://arxiv.org/abs/gr-qc/9405047
This also seems to be related, but more recent; it gauges the Bargmann group and also works on its relativistic version (that I've essentially been describing above), which it called the "extended Poincaré group"
A pedagogical review of gravity as a gauge theory
https://arxiv.org/abs/2104.02627
Edit: I also just found this reference Review on Non-Relativistic Gravity that seems to be directly addressing and resolving the issues I raised. It's got the 1984 Duval, Burdet, Künzle and Perrin reference in it, as well as the De Pietri, Lusanna and Pauri, so it's working off of them, too. Take a look at what they have going on in sections 5 and 6, and compare their version of the non-relativistic Schwarzschild in section 7.2 with what I put up. They're using torsion to get around some of the limitations of non-relativistic gravity.
