Is Gauss law for gravity local? in Newtonian gravity, the gravitational field obeys the equation
$$\nabla^2 \phi = 4 \pi G \rho$$
David Tong in his notes on general relativity claims that this equation works well when $\rho$ is not a function of time, but when $\rho$ depends on time, then a change in $\rho$ would result in a change in $\phi$ throughout all of space. But this equation is a local equation, in the sense that I'm calculating the Laplacian of $\phi$ at the point $\vec{r}$ at the time $t$, and this is proportional to the density at the position $\vec{r}$ at the time $t$. So  how does a change in $\rho$ at one position affect the field at another position instantaneously according to this equation?
 A: Recall, that given a point with mass $M$, its gravitational potential depending on the distance $r$ to the point is given by:
$$\phi(r)=-G\frac{M}{r}.$$
Recall as well, that given a continuum with density $\rho$, its gravitational potential at $\vec{r}$ is given by integrating over all the contributions:
$$\Phi(t,\vec{r})
=-G\int\mathrm d^3\vec{r}'\frac{\rho(t,\vec{r}')}{|\vec{r}-\vec{r}'|}.$$
The problem here, is that when you move from the point $\vec{r}$ to the point $\vec{r}'$ to gather information about its contribution to the gravitational field (at $\vec{r}$!), you should also move into the past as the information took a time of $|\vec{r}-\vec{r}'|/c$ to travel from $\vec{r}'$ to $\vec{r}$. In the equation above, which can be transfered into your differential equation by applying the Laplace operator and using its fundamental solution:
$$-\Delta_{\vec{r}}\frac{1}{|\vec{r}-\vec{r}'|}=4\pi\delta^3(\vec{r}-\vec{r}')$$
does not consider that. Changing $\rho(t,\vec{r}')$ for any $\vec{r}'$ does affect $\Phi(t,\vec{r})$, even if $\vec{r}$ and $\vec{r}'$ are many light-years apart. The only thing that matters for this to work is the same time of $t$ on both sides. In the integral, every contribution from everywhere in space is regarded at this current time $t$, implying that the information travels with infinite speed. A better approach is therefore the following generalized integral using the speed of light $c$ with which the information about the gravitational field travels:
$$\Phi(t,\vec{r})
=G\int\mathrm d^3\vec{r}'\frac{\rho\left(t-\frac{|\vec{r}-\vec{r}'|}{c},\vec{r}'\right)}{|\vec{r}-\vec{r}'|},$$
which can be transfered into a differential equation (the wave equation) by applying the D'Alembert operator (wave operator) and using its fundamental solution:
$$\square_{t,\vec{r}}\frac{\delta^1\left(t-t'-\frac{|\vec{r}-\vec{r}'|}{c}\right)}{|\vec{r}-\vec{r}'|}
=4\pi\delta^1(t-t')\delta^3(\vec{r}-\vec{r}').$$
The following wave equation now includes the effect of retardation (it is really called retarded gravity) and falls back to the Laplace equation either if $c\rightarrow\infty$ or if $\rho$ (and therefore $\Phi$) are independent of time like you mentioned, that David Tong claimed:
$$\square\Phi
=\frac{1}{c^2}\frac{\partial^2}{\partial t^2}\Phi-\Delta\Phi
=-4\pi G\rho.$$
By the way, changing the sign in front of $|\vec{r}-\vec{r}'|/c$ from $-$ to $+$ changes the retarded potential into an advanced potential and instead of information from the past, information from the future enters your integral, which is why advanced potentials are rarely considered in physics. But they are: Retarded potentials describe the emission while advanced potentials describe the absorption of fields.
A: In some cases, the value of the field in the neighbourhood of some point $\mathbf r_0$ is such that $\nabla^2\phi = 4\pi G \rho(\mathbf r_0)$
It happens for example outside a massive sphere, where $\rho = 0$, or inside a massive sphere with constant density. But it is not the case for a massive sphere with non constant density. For a point inside a hole from pole to pole across a sphere, with hole diameter much smaller than the sphere radius, $\nabla^2\phi \neq 0$, in spite of $\rho = 0$ locally.
In general, it is necessary to apply the Green function for this differential equation, which is an integral over all of space.
$$\phi(\mathbf{r}) = - \iiint  \frac{\rho(\mathbf{r}')}{ |\mathbf{r} - \mathbf{r}'|}\, \mathrm{d}^3\! r'$$
