I'm a student, trying to re-derive a result found in a paper by calculating the following in spherical coordinates:

$$\mathbf{I}+ \frac{\nabla\nabla}{constant}$$

where $\mathbf{I}$ is a $3x3$ identity matrix.

The paper that I've seen writes the result as

\begin{bmatrix} 0 & 0 & 0\\ cos(\Theta)cos(\Phi) & cos(\Theta)sin(\Phi) & -sin(\Theta)\\ -sin(\Phi) & cos(\Theta) & 0 \end{bmatrix}

How do they get that? What is the gradient of a gradient in spherical coordinates? Is it a Hessian?

The $\nabla \nabla$ here is not a Laplacian (divergence of gradient of one or several scalars) or a Hessian (second derivatives of a scalar), it is the gradient of the divergence. That is why it has matrix form: it takes a vector and outputs a vector. (Taking the divergence of a vector gives a scalar, another gradient yields a vector again).