$$ \textbf{F}_{grav} = m \textbf{g} = - \left( \frac{GmM}{r^2} \right) \textbf{e}_r = -m \boldsymbol{\nabla} \Phi (\textbf{x}) $$

where $ \textbf{e}_r $ is the unit vector pointing from the first particle to the other and $ \Phi (\textbf{x}) $ is the gravitational potential from the first particle. Solve for the acceleration in terms of the gravitational potential and you get $ \textbf{a} \equiv \textbf{g} \equiv - \boldsymbol{\nabla} \Phi $ . You can rewrite this using index notation to make it more compact. It should then look like this

$$ \frac{d^2 x^i}{dt^2} = - \eta^{ij} \frac{\partial \Phi}{\partial x^j} = - \eta^{ij} [\partial_j \Phi]_\textbf{x} $$

Here $ x^i $ is the position of the particle in our frame of reference and $ \eta^{ij} $ represents the spatial components of flat spacetime. This doesn't change the equation but only helps us in balancing out our indices. The equation above tells us how one of our particles falls, but now we need to write a similar equation for the second particle. However, this second particle is at some distance $ \textbf{n} (t) $ away from the first particle; we call this a separation vector. This means that the position of the second particle is at $ x^i(t) + n^i(t) $. With this information we can now write the equation of motion for the second particle, which is

$$ \frac{d^2 (x^i + n^i)}{dt^2} = - \eta^{ij} [\partial_j \Phi]_\textbf{x + n} $$

Now we must expand the right side in a Taylor series in order to get an estimate of the derivative. This is what the expansion around the reference looks like

$$ [\partial_j \Phi]_\textbf{x + n} = [\partial_j \Phi]_\textbf{x} + n^k \left( \frac{\partial}{\partial x^k} [\partial_j \Phi] \right)_\textbf{x} + \cdots $$

This works assuming that the quantity of $ \textbf{n} $ is small, so we ignore higher order terms. Now, we are able to substitute this expansion into our equation of motion for the second particle, and then subtract the equation of motion of the first particle from that of the second particle.

1) Plug in the Taylor series expansion into the equation of our second particle as follows

$$ \frac{d^2 (x^i + n^i)}{dt^2} = - \eta^{ij} [\partial_j \Phi]_\textbf{x} - \eta^{ij} n^k \left( \frac{\partial}{\partial x^k} [\partial_j \Phi] \right)_\textbf{x} $$

where

$$ \frac{d^2 (x^i + n^i)}{dt^2} = \frac{d^2 x^i}{dt^2} + \frac{d^2 n^i}{dt^2} $$

2) Subtract the equation of the first particle from that of the second like this

$$ \frac{d^2 x^i}{dt^2} + \frac{d^2 n^i}{dt^2} - \frac{d^2 x^i}{dt^2} = - \eta^{ij} [\partial_j \Phi]_\textbf{x} - \eta^{ij} n^k [\partial_k \partial_j \Phi]_\textbf{x} + \eta^{ij} [\partial_j \Phi]_\textbf{x} $$

3) Simplify by canceling out terms to get

$$ \frac{d^2 n^i}{dt^2} \approx - \eta^{ij} [\partial_k \partial_j \Phi]_\textbf{x} n^k \equiv - \eta^{ij} \left( \frac{\partial^2 \Phi}{\partial^k \partial^j} \right) n^k $$

This is the

**Newtonian Deviation Equation.**By knowing the separation of the two particles at one time and assuming the separation vector remains small, we can determine the deviation of the particles at any given time.