The formula looks a bit unusual but note that
\begin{align}
\frac1n\sum_{i=1}^n(S_i-\bar S)(S_i-\bar S)^T
&=\frac1n\sum_{i=1}^n (S_iS_i^T-2\bar S S_i^T+\bar S \bar S^T)
\\&=\frac1n\sum_{i=1}^n S_iS_i^T-2\bar S\frac1n\sum_{i=1}^n S_i^T+\frac1n\bar S\bar S^T
\\&=\frac1n\sum_{i=1}^n S_iS_i^T-\frac1{n^2}\sum_{i=1}^nS_i\sum_{i=1}^nS_i^T.
\end{align}
This can be seen as a kind of estimator of $\operatorname{Var}S_i$. If the observations (and the contributions $S_i$ to the score vector) are independent, $n$ times this is a kind of estimator of the expected (and I suppose the observed) Fisher information
$$
I(\theta)=\operatorname{Var}S=\operatorname{Var}(\sum_{i=1}^n S_i)=n\operatorname{Var}S_i.
$$
But if the likelihood is available, the observed Fisher information is directly observable so it is a bit strange that the authors try to only estimate it in this way.
Note also that the last term in their expression is not needed as $S=\sum S_i=0$ at the MLE of $\theta$.
This estimator (without the last term) is known as the empirical Fisher information in the machine learning literature. Kunstner et. al. 2019 argue "that the conditions under which the empirical Fisher approaches the Fisher (and the Hessian) are unlikely to be met in practice, and that, even on simple optimization problems, the pathologies of the empirical Fisher can have undesirable effects."
