  We consider Steklov eigenvalues on nearly spherical and nearly annular
domains in $d$ dimensions. By using the Green-Beltrami identity for spherical
harmonic functions, the derivatives of Steklov eigenvalues with respect to the
domain perturbation parameter can be determined by the eigenvalues of a matrix
involving the integral of the product of three spherical harmonic functions. By
using the addition theorem for spherical harmonic functions, we determine
conditions when the trace of this matrix becomes zero. These conditions can
then be used to determine when spherical and annular regions are critical
points while we optimize Steklov eigenvalues subject to a volume constraint. In
addition, we develop numerical approaches based on particular solutions and
show that numerical results in two and three dimensions are in agreement with
our analytic results.
