The Wikipedia page on the Lamb shift includes the following first steps:
$$\Delta V = V\bigl(\vec{r}+\delta \vec{r}\bigr)-V(\vec{r})=\delta \vec{r} \cdot \nabla V (\vec{r}) + \frac{1}{2} \bigl(\delta \vec{r} \cdot \nabla\bigr)^2V(\vec{r})+\cdots\tag{1}$$
Because the fluctuations are isotropic:
$$\begin{align} \langle \delta \vec{r} \rangle _\text{vac} &=0 \tag{2.1}\\ \langle (\delta \vec{r} \cdot \nabla )^2 \rangle _\text{vac} &= \frac{1}{3} \langle (\delta \vec{r})^2\rangle _\text{vac} \nabla ^2\tag{2.2} \end{align}$$
Then
$$\langle \Delta V\rangle =\frac{1}{6} \langle (\delta \vec{r})^2\rangle _\text{vac}\left\langle \nabla ^2\left(\frac{-e^2}{4\pi \epsilon _0r}\right)\right\rangle _{at}\tag{3}$$
I understand that step (1) is a Taylor expansion. Step (2.1) also makes intuitive sense, isotropy means that it is the same in every direction so it's no surprise that the expectation of $\delta \vec{r}$ is zero.
I can also understand step (3), which seems to be substituting results from step (2) into step (1), while ignoring higher order terms.
However, I have no clue where step (2.2) comes from. I attempted naively expanding the dot product as you would with $(a \cdot b)^2$, but I don't know if this is allowed.
In the same vein, notation wise, is $\delta \vec{r} \cdot \nabla$ the same as $\nabla \cdot \delta \vec{r}$? I know that the scalar product is commutative, but then $\nabla \cdot (\text{something})$ returns the divergence.
TL;DR: Don't know where step 2 comes from, and getting confused by vector calculus.