From pg. 35 of Classical Electrodynamics 3rd edition, Jackson,
$$\begin{aligned} \nabla^{2} \Phi_{a}(\mathbf{x}) &=-\frac{1}{4 \pi \epsilon_{0}} \int \rho\left(\mathbf{x}^{\prime}\right)\left[\frac{3 a^{2}}{\left(r^{2}+a^{2}\right)^{5 / 2}}\right] d^{3} x^{\prime} \end{aligned}$$
"Choose R such that $\rho(\mathbf{x'})$ changes little over the interior of the sphere... With a Taylor series expansion of the well-behaved $\rho (\mathbf{x'})$ around $\mathbf{x'} = \mathbf{x}$ one finds ..."
\begin{align} \nabla^{2} \Phi_{a}(\mathbf{x}) &=-\frac{1}{\epsilon_{0}} \int_{0}^{R} \frac{3 a^{2}}{\left(r^{2}+a^{2}\right)^{5 / 2}}\left[\rho(\mathbf{x})+\frac{r^{2}}{6} \nabla^{2} \rho+\cdots\right] r^{2} d r+O\left(a^{2}\right), \end{align}
where $r = |\mathbf{x'} -\mathbf{x}|$.
Could someone explain how to derive this Taylor series for a function of a vector? I've never seen this before and am at a loss.