The other day, I asked myself the question: How would you write down Newtonian gravity as a classical field theory?
The goal, of course, is to write down a Lagrangian density that yields, as the corresponding Euler-Lagrange equation, Poisson's equation for gravity:
\begin{align}
\nabla^2\phi=4\pi G\rho,
\end{align}
where $\phi$ is the gravitational potential field and $\rho$ is the mass density.
The seemingly obvious choice for a Lagrangian is given by
\begin{align}
{\cal L}=\frac{1}{8\pi G}(\nabla\phi)^2+\phi\rho,
\end{align}
as variation of this Lagrangian with respect to $\phi$ indeed yields Poisson's equation. This Lagrangian is even presented in some Wikipedia articles on the topic. Unfortunately, it is not the complete picture: since $\rho$ is also a dynamical field, it must also be varied, and variations with respect to $\rho$ yield the rather nonsensical equation,
\begin{align}
\phi=0.
\end{align}
Oops.
Edit (Nov 11, 2015): Consider what follows below a whimsical flight of fancy that is really misguided: by trying to construct a "matter" part for the gravitational Lagrangian, I made the mistake of assuming that I in fact have a theory of matter. The proper thing to do would have been to leave the matter part unspecified: \({\cal L}=\frac{1}{8\pi G}(\nabla\phi)^2+{\cal L}_M\), prescribing only that its variation with respect to \(\phi\) must yield \(\rho\): \(\delta{\cal L}_M/\delta\phi=\rho\). Poisson's equation automatically follows.
How can this be fixed? Obviously, the second equation should capture the dynamics of how matter responds to gravity, just as the first equation told us what gravitational field is induced by matter. Now we know that matter actually responds to the gradient of the gravitational field. In fact, what we seek (for collisionless dust in the absence of pressure, viscosity, and other nasty things that would require the full machinery of fluid dynamics) is the equation
\begin{align}
\vec{a}=-\nabla\phi,
\end{align}
i.e., the standard acceleration law in the presence of gravity for the acceleration $\vec{a}$.
This requirement immediately suggests, however, a Lagrangian density in the form
\begin{align}
{\cal L}=\frac{\rho_0}{8\pi G\rho}(\nabla\phi)^2+\phi\rho_0+\frac{\rho_0}{8\pi G}\int\frac{a^2}{\rho^2}\dot{\rho}dt,
\end{align}
where $t$ is time and the overdot is differentiation with respect to time.
Now a few paragraphs above, I insisted that $\rho$ must be treated as a dynamical field and varied accordingly. Why is the same thing not true for $a$? Because it is not a new degree of freedom: it is just the acceleration associated with $\rho$. In other words, as we derive the Euler-Lagrange equation from this Lagrangian, we treat $a$ as though it was a function of $\rho$ (a multivalued function, to be precise, or rather, a multivalued relation or mapping, but the point is, it is not an independent degree of freedom.)
Variation of this Lagrangian with respect to $\rho$ yields the equation,
\begin{align}
-\frac{\rho_0}{8\pi G}\frac{(\nabla\phi)^2}{\rho^2}+\frac{\rho_0}{8\pi G}\frac{a^2}{\rho^2}=0,
\end{align}
or
\begin{align}
|\vec{a}|=|\nabla\phi|,
\end{align}
which is almost the equation we sought. It is consistent with $\vec{a}=-\nabla\phi$ though other directions for $\vec{a}$ are not excluded so long as its magnitude remains the same. In principle, we could ensure that these quantities remain parallel by introducing it as a constraint through a Lagrange-multiplier term, such as $\lambda(\vec{a}\times\nabla\phi)^2$, which does not alter the preceding two field equations but does enforce the parallel constraint, but this goes beyond the point I am trying to get to here.
Which is that there is a reason why I was interested in this approach. I was studying non-local Lagrangians. Non-local Lagrangians tend to have terms that contain integrals that are taken over all of space. The Green's function solution to Poisson's equation is just such an integral. What if we were to use it in place of $\phi$ in the above Lagrangian?
The solution to Poisson's equation is given by
\begin{align}
\phi(\vec{r})=-G\int\frac{\rho(\vec{x}')}{|\vec{x}-\vec{x}'|}d^3\vec{x}',
\end{align}
so the Lagrangian then reads,
\begin{align}
{\cal L}=\frac{G\rho_0}{8\pi\rho}\left[\nabla\int\frac{\rho(\vec{x}')}{|\vec{x}-\vec{x}'|}d^3\vec{x}'\right]^2-G\rho_0\int\frac{\rho(\vec{x}')}{|\vec{x}-\vec{x}'|}d^3\vec{x}'+\frac{\rho_0}{8\pi G}\int\frac{a^2}{\rho^2}\dot{\rho}dt.
\end{align}
The second term in this Lagrangian density is obviously a surface term and thus, can be omitted, leaving us with
\begin{align}
{\cal L}=\frac{G\rho_0}{8\pi\rho}\left[\nabla\int\frac{\rho(\vec{x}')}{|\vec{x}-\vec{x}'|}d^3\vec{x}'\right]^2+\frac{\rho_0}{8\pi G}\int\frac{a^2}{\rho^2}\dot{\rho}dt,
\end{align}
and variation with respect to $\rho$ yields, after getting rid of further surface terms, the equation,
\begin{align}
|\vec{a}|=\left|G\nabla\int\frac{\rho(\vec{x}')}{|\vec{x}-\vec{x}'|}d^3\vec{x}'\right|,
\end{align}
which is the non-local, "action-at-a-distance" version of Newton's gravity.
This now is a non-local Lagrangian of only one dynamical field. What this suggests is that it is possible to replace a dynamical field in a Lagrangian with a non-local term. Conversely, I guess it is also possible that formally non-local terms in a Lagrangian can be made strictly local, if the non-local term is the solution of a field equation, by introducing a new field variable.
