Consider an instrument value $f(S_0^1, \ldots, S_0^n)$ that depends on $n$ spot levels. Let $$\overrightarrow{S}_0=[S_0^1, \ldots, S_0^n]^T$$ be an $n$-dimensional vector representing the spot levels. We can approximate the cross gamma
\begin{align*}
\frac{\partial^2 f\big(\overrightarrow{S}_0\big)}{\partial S_0^i \partial S_0^j}
\end{align*}
by a finite difference scheme of the form
\begin{align*}
&\frac{1}{4 \varepsilon^2 S_0^i S_0^j}\Big[f\big(\overrightarrow{S}_0 + \varepsilon S_0^i\overrightarrow{1}_i + \varepsilon S_0^j\overrightarrow{1}_j\big) - f\big(\overrightarrow{S}_0 + \varepsilon S_0^i\overrightarrow{1}_i - \varepsilon S_0^j\overrightarrow{1}_j\big)\\
& \qquad\qquad - f\big(\overrightarrow{S}_0 - \varepsilon S_0^i\overrightarrow{1}_i + \varepsilon S_0^j\overrightarrow{1}_j\big) + f\big(\overrightarrow{S}_0 - \varepsilon S_0^i\overrightarrow{1}_i - \varepsilon S_0^j\overrightarrow{1}_j\big)\Big],
\end{align*}
where $\overrightarrow{1}_i$ is an $n$-dimensional vector with $1$ in the $i^{th}$ element and zeros elsewhere. Here, $\varepsilon$ is a small perturbation, which is usually set to $0.01$.
To justify the above, we use Taylor expansion. For example,
\begin{align*}
f\big(\overrightarrow{S}_0 + \Delta_i\overrightarrow{1}_i + \Delta_j\overrightarrow{1}_j\big)&\approx f\big(\overrightarrow{S}_0\big)+\Big(\Delta_i\frac{\partial }{\partial S_0^i}+ \Delta_j\frac{\partial }{\partial S_0^j}\Big)f\big(\overrightarrow{S}_0\big)\\
& \ \ + \frac{1}{2}\Big(\Delta_i\frac{\partial }{\partial S_0^i}+ \Delta_j\frac{\partial }{\partial S_0^j}\Big)^2f\big(\overrightarrow{S}_0\big) + o(\Delta_i \Delta_j ).
\end{align*}