Currently working my way through Gravity Newtonian, Post-Newtonian, Relativistic - Poisson and Will (2014) and I've stumbled upon a potential issue with my memory of vector calc in index notation.
Context
I've got an equation which takes the following form
$$ \rho ( \partial_t \mathbf{v} + \mathbf{v} \cdot \nabla \mathbf{v}) = \rho \nabla U - \nabla p, \tag{1} \label{eq:maineq} $$
where mass density, pressure, velocity vector field and Newtonian gravitational potential for a fluid element are given by $\rho, p, \mathbf{v}$ and $U$ respectively. Usual notation is used i.e. $\partial_t\equiv\partial/\partial t$ etc. etc.
Now, in terms of a fluid element comprised of a continuous matter distribution, the Newtonian gravitational potential satisfies Poisson's equation, namely
$$ \nabla^2U = -4\pi G \rho. \tag{2} $$
We can express the first term on the RHS of Eq. (\ref{eq:maineq}) as the following
$$ \rho \partial_j U = -\frac{1}{4 \pi G} \nabla^2 U \partial_j U. $$
where $\partial_j U \equiv \nabla U$. The steps that follow is where the issue lies.
Question
We have $ -\frac{1}{4 \pi G} \nabla^2 U \partial_j U$ which can be expressed as the following
$$ -\frac{1}{4 \pi G} (\partial_k\partial_k U) \partial_j U. \tag{3} \label{eq:laplacegrad}$$
In light of the comments below can someone explain the steps that Poisson and Will have taken to arrive at the following:
$$\rho \partial_j U = -\frac{1}{4 \pi G} \partial_k \left( \partial_j U \partial_k U - \frac{1}{2} \delta^{jk} \partial_i U \partial_i U \right). \tag{4} \label{eq:derive}$$
With regards the rules of the site. To see my (incorrect) attempt, it can be seen in a previously edited version of the question.
See below for answer.
Answer ( thanks to the comments above)
To summarise the comments above and the steps necessary for the correct derivation, we may begin by noticing
$$(\nabla^2 U) \partial_j U \equiv (\partial_k \partial_k U) \partial_j U.$$
Apply the product rule to $\partial_k (\partial_k U \partial_j U)$ giving
$$ \partial_k (\partial_k U \partial_j U)= (\partial_k\partial_k U)\partial_j U + \partial_k U (\partial_k \partial_j U), $$
where we may deduce
$$ (\partial_k \partial_k U) \partial_j U = \partial_k (\partial_k U \partial_j U) - \partial_k U (\partial_k \partial_j U). $$
Next, notice the relation
$$ \partial_k U (\partial_k \partial_j U) \equiv \frac{1}{2} \partial_j |\nabla U|^2. $$
Which can be seen by performing the chain rule
$$ \frac{1}{2} \partial_j |\nabla U|^2 = \frac{1}{2} \cdot 2 \partial_k U (\partial_k \partial_j U). $$
Finally, we can express $\partial_j \equiv \partial_k \delta^{k}_{j}$ and then piecing together allows one to successfully derive the relation given in the question above.