The solution is straightforward in Cartesian coordinates and impossibly messy in spherical coordinates. I will describe it for a cone that fits within a hemisphere; larger cones can be treated in a similar manner.
Choose units in which the sphere has unit radius. After a suitable rotation by some orthogonal matrix $\mathbb{Q}$, the cone will be vertical and projects onto a circle of radius $\rho$ in the xy plane. In cylindrical coordinates $(r, \theta, z)$, the points will have a uniform distribution in $\theta$ between $0$ and $2\pi$ (by cylindrical symmetry), a uniform distribution in $z$ between $\sqrt{1-\rho^2}$ and $1$, and necessarily $r = \sqrt{1-z^2}$.
Because $z$ has a uniform distribution, its variance is
$$(1 - (1-\rho^2))/12 = \rho^2/12.$$
The distribution of $r$ is deduced from that of $z$ by taking the Jacobian; its pdf is
$$\frac{r}{\sqrt{1-r^2} (1 - \sqrt{1-\rho^2})}, \quad 0 \le r \le \sqrt{1-\rho^2}.$$
By axial symmetry the expectations of $x$ and $y$ are both zero, whence their variances are the expectations of their squares. By the reflection symmetry $x\to y$, $y\to x$, these expectations are equal. The sum of those expectations is the expectation of $x^2+y^2 = r^2$, which can be computed as
$$\int_0^{\rho } \frac{ r^2 \, rdr}{\left(1-\sqrt{1-\rho ^2}\right) \sqrt{1-r^2}} = \frac{2- \left(2+\rho ^2\right)\sqrt{1-\rho ^2} }{3(1-\sqrt{1-\rho ^2})}.$$
Half of this quantity therefore is the common variance of $x$ and $y$. (As a quick check, the limiting values of the $x$ and $y$ covariances at $\rho\to 1$ are $1/3$. This is correct, because the variance of $x$ for the upper hemisphere equals the variance of $x$ for the sphere, which clearly is $1/3$.)
The symmetries $x\to -x$ and $y\to -y$ immediately imply the covariances among $x$, $y$, and $z$ are all $0$. We have thereby obtained the full covariance matrix $\mathbb{D}$. Applying the original rotation gives the original covariance matrix as $\mathbb{Q'DQ}$.
