  We consider a topological insulator (TI) of spherical geometry and
numerically investigate the influence of disorder on the density of surface
states. To the clean Hamiltonian we add a surface disorder potential of the
most general hermitian form, $V = V^0(\theta,\phi) + {\bf V}(\theta,\phi) \cdot
{\bf \sigma}$. We expand these four disorder functions in spherical harmonics
and draw the expansion coefficients randomly from a four-dimensional, zero-mean
gaussian distribution. Different strengths and classes of disorder are realized
by specifying the $4 \times 4$ covariance matrix. For each instantiation of the
disorder, we solve for the energy spectrum via exact diagonalization. Then we
compute the disorder-averaged density of states, $\rho(E)$, by averaging over
200,000 different instantiations. Disorder broadens the Landau-level delta
functions of the clean density of states into peaks that decay and merge
together. If the spin-dependent term is dominant, these peaks split due to the
breaking of the degeneracy between time-reversed partner states. Increasing
disorder strength pushes states closer and closer to zero energy (the Dirac
point), resulting in a low-energy density of states that becomes nonzero for
sufficient disorder, typically approaching an energy-independent saturation
value, for most classes of disorder. But for purely spin-dependent disorder
with ${\bf V}$ either entirely out-of-surface or entirely in-surface, we
identify intriguing disorder-induced features in the vicinity of the Dirac
point. In the out-of-surface case, a new peak emerges at zero energy. In the
in-surface case, we see a symmetry-protected zero at zero energy, with
$\rho(E)$ increasing linearly toward nonzero-energy peaks. These striking
features are explained in terms of the breaking (or not) of two chiral
symmetries of the clean Hamiltonian.
