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.
