I'm trying to derive the second functional derivative of the Weizsäcker kinetic energy functional in 3D given by $T_\text{vW}[n] = \int d^3\mathbf{r} \frac{1}{8} \frac{|\nabla n(\mathbf{r})|^2}{n(\mathbf{r})}$, with boundary conditions where $n(\mathbf{r})\to0$ at infinity.
My exact problem, using the the method described below, is fundamentally a vector-calculus manipulation problem if one would like to skip the "Preamble/Context" directly to "My Problem". Otherwise, it's a functional-calculus problem if one would rather suggest alternative methods that might be more straightforward to get the second functional derivative of such expressions involving 3D gradients.
Preamble/Context
My strategy is to use the following relations:
$\frac{dT[n+\epsilon\phi]}{d\epsilon} \Bigr\rvert_{\epsilon = 0} = \int d^3\mathbf{r} \frac{\delta T[n]}{\delta n(\mathbf{r})} \phi(\mathbf{r})$
$\frac{d^2T[n+\epsilon\phi]}{d\epsilon^2} \Bigr\rvert_{\epsilon = 0} = \int d^3\mathbf{r}~d^3\mathbf{r}' \frac{\delta^2 T[n]}{\delta n(\mathbf{r})\delta n(\mathbf{r}')} \phi(\mathbf{r}) \phi(\mathbf{r}')$
i.e. plug the functional into the LHS and manipulate until it looks like the RHS and extract the functional derivatives by inspection. For example, for the first functional derivative,
$\begin{equation} \begin{aligned} \frac{dT_\text{vW}[n+\epsilon\phi]}{d\epsilon} \Bigr\rvert_{\epsilon = 0} &= \frac{1}{8} \left[- \int d^3\mathbf{r} \frac{|\nabla n(\mathbf{r}) +\epsilon\nabla \phi(\mathbf{r})|^2}{(n(\mathbf{r}) +\epsilon\phi(\mathbf{r}))^2} \phi(\mathbf{r}) + \int d^3\mathbf{r} \frac{2\nabla n(\mathbf{r}) \cdot\nabla\phi(\mathbf{r})+2\epsilon|\nabla\phi(\mathbf{r})|^2}{n(\mathbf{r})+\epsilon \phi(\mathbf{r})} \right]\Bigr\rvert_{\epsilon = 0} \\ &= \frac{1}{8} \left[- \int d^3\mathbf{r} \frac{|\nabla n(\mathbf{r}) |^2}{n^2(\mathbf{r})} \phi(\mathbf{r}) + 2 \int d^3\mathbf{r} \frac{\nabla n(\mathbf{r})}{n(\mathbf{r})} \cdot\nabla\phi(\mathbf{r}) \right] \\ &= \frac{1}{8} \left[- \int d^3\mathbf{r} \frac{|\nabla n(\mathbf{r}) |^2}{n^2(\mathbf{r})} \phi(\mathbf{r}) + 2 \int \phi(\mathbf{r}) \frac{\nabla n(\mathbf{r})}{n(\mathbf{r})} \cdot d\mathbf{S} - 2\int d^3\mathbf{r} \nabla \cdot \left(\frac{\nabla n(\mathbf{r})}{n(\mathbf{r})}\right) \phi(\mathbf{r}) \right] \\ &= \int d^3\mathbf{r} \left[\frac{1}{8} \frac{|\nabla n(\mathbf{r}) |^2}{n^2(\mathbf{r})} - \frac{1}{4} \frac{\nabla^2 n(\mathbf{r})}{n(\mathbf{r})} \right] \phi(\mathbf{r}) \end{aligned} \end{equation}$
where I have assumed that $\phi(\mathbf{r})$ must obey the same boundary conditions (going to zero at infinity) and have used the boundary conditions to eliminate surface integral terms. By comparing with (1), we see that $\frac{\delta T[n]}{\delta n(\mathbf{r})} = \frac{1}{8} \frac{|\nabla n(\mathbf{r}) |^2}{n^2(\mathbf{r})} - \frac{1}{4} \frac{\nabla^2 n(\mathbf{r})}{n(\mathbf{r})}$
(The Wikipedia article on functional derivatives uses a different method but gets the same result)
For the second functional derivative, the same process was repeated, but with a second derivative wrt $\epsilon$ this time. If I'm not mistaken, plugging in the functional trivially results in
$\frac{d^2T[n+\epsilon\phi]}{d\epsilon^2} \Bigr\rvert_{\epsilon = 0} = \frac{1}{4} \int d^3\mathbf{r} \frac{|\nabla n(\mathbf{r})|^2 \phi^2(\mathbf{r})}{n^3(\mathbf{r})} + \frac{|\nabla \phi(\mathbf{r})|^2}{n(\mathbf{r})} - 2 \frac{\nabla n(\mathbf{r})\cdot \nabla \phi(\mathbf{r})\phi}{n^2(\mathbf{r})}$
Terms with $\phi^2(\mathbf{r})$ can then be treated easily with
$\phi^2(\mathbf{r}) = \int d^3\mathbf{r}' \phi(\mathbf{r})\phi(\mathbf{r}')\delta(\mathbf{r}-\mathbf{r}')$
to match the form required. Hence, the goal now is to rewrite the $\nabla \phi$ terms in terms of $\phi$ only. After fiddling around with some vector calculus identities, I managed to express the third term as
$ - 2 \int d^3\mathbf{r} \frac{\nabla n(\mathbf{r})\cdot \nabla \phi(\mathbf{r})\phi}{n^2(\mathbf{r})} = \int d^3\mathbf{r} \left(\frac{\nabla^2 n(\mathbf{r})}{n^2(\mathbf{r})} - 2 \frac{|\nabla n(\mathbf{r})|^2}{n^3(\mathbf{r})}\right)\phi^2(\mathbf{r})$
My Problem
I can't figure out how to rewrite the second term $\int d^3\mathbf{r} \frac{|\nabla \phi(\mathbf{r})|^2}{n(\mathbf{r})}$ solely in terms of $\phi(\mathbf{r})$ and not any of its higher derivatives in order to match the required form in (2). I also couldn't find the form of the desired second functional derivative online, so I'm stuck not even knowing the form I'm supposed to force it into!
Any help is much appreciated, thank you!