Investigate maxima of Gaussian integral over sphere. Let $\alpha>0$ be a positive parameter and consider the function
$$f(x) = \int_{\mathbb S^{n-1}} e^{-\alpha \left\lVert x-y \right\rVert^2} dS(y)$$
for $x \in \mathbb R^n.$ So, since this was asked, although we integrate over the unit sphere, the function "lives" on $\mathbb R^n.$
This function is clearly rotationally symmetric. 
I would like to show that the global maxima are attained at one single radius $r$, only. 
The rotational symmetry implies that we can consider it as a one-dimensional function by choosing $x=(x_1,0....,0)$, this way the exponent simplifies to $e^{-\alpha \left\lVert x-y \right\rVert^2}=e^{-\alpha (x_1-y_1)^2+1-y_1^2}.$
If anything is unclear about this question, then please let me know. I am happy to hear about any ideas how to approach this problem.
EDIT: Thanks to some interesting comments below, one can say that the global maximum is always attained at some radius $r \in [0,1]$ where for small $\alpha$ it seems to be attained close to zero and for large $\alpha$ it is attained closer to one. 
The question remains however why is there only one radius at which the global maximum is attained? -In fact as George Lowther points out in the comments, for $\alpha \le n/2$ the unique maximum is attained at $r=0$ which leaves the case $\alpha >n/2$ when this does not hold true.
 A: This is certainly not a full solution.
First we would like to get an expression for $f(x)$.  (Below I write $f(x)$ or $f_n(x)$ when I should probably should have written $f((x_1, 0,0,\ldots, 0)$.)
I am going to use the function
$$
s(n,r) := \frac{2 \pi ^{\frac{n+1}{2}}}{\Gamma \left(\frac{n+1}{2}\right)} r^n,
$$
the surface area of an n-dimensional sphere. 
We can use $s$ to integrate $f$,
$$
f_n(x) = \int_{\theta=0}^\pi s (n-1 ,\sin(\theta)) \exp\left(-a \left((x-\cos (\theta))^2+\sin ^2(\theta)\right)\right)d\theta.
$$
According to Mathematica, $\;f_2(x) = \frac{\pi  e^{-a (x+1)^2} \left(e^{4 a x}-1\right)}{a x}$.
And $f_4(x) = \frac{\pi ^2 e^{-a (x+1)^2} \left(2 a x+e^{4 a x} (2 a x-1)+1\right)}{2 a^3 x^3}$.
For odd values of $n$, Mathematica gives Hypergeometric functions. 
$$f_3(x) = 2 \pi ^2 e^{-a \left(x^2+1\right)} \, _0\tilde{F}_1\left(;2;a^2 x^2\right)$$
In the Mathematica language, $f_3(x)$ is 
"(2*Pi^2*Hypergeometric0F1Regularized[2, a^2*x^2])/E^(a*(1 + x^2))".
If you try to find the maximum of $f$, it seems that the maximum occurs at $x=0$ when $a<2$.  One way to think about this is releasing a million drunkards randomly on the $S^{n-1}$ sphere and having them walk randomly.  Initially, the most drunkards per unit volume occurs on the $S^{n-1}$ sphere, but as time goes on the radius of maximum drunkards moves toward the center.  
Another way to look at it is as a convolution of the sphere with a Gaussian with $\sigma = \sqrt{1/(2 a)}$.  When sigma is large, $f_n(x)$ is approximately Gaussian.
I was not able to find closed form solutions for $\mathrm{argmax}_x\, f_n(x)$.  If we take the derivative when $n$ is even, set it equal to zero, and simplify, we get expressions of the form 
$$p_1(x,a) + \cosh(2 a x)p_2(x,a) +\sinh(2 a x)p_2(x,a) =0$$
where the $p_i$ are polynomials of degree of low degree with "nice" coefficients.  
For example, for $n=10$, Mathematica says that the derivative of $f_{10}$ is zero if 
$$32 a^5 x^6 \sinh (2 a x)-32 a^5 x^5 \cosh (2 a x)-160 a^4 x^5 \cosh (2 a x)+240 a^4
   x^4 \sinh (2 a x)+360 a^3 x^4 \sinh (2 a x)-840 a^3 x^3 \cosh (2 a x)-420 a^2 x^3 \cosh (2 a
   x)+1680 a^2 x^2 \sinh (2 a x)+210 a x^2 \sinh (2 a x)+945 \sinh (2 a x)-1890 a x \cosh (2 a
   x)=0.$$
Hope that helps.
A: If $\mathcal S$ is the unit $(n - 1)$-sphere around the origin and $\lVert x \rVert = r$, then choosing hyperspherical coordinates s.t. $x_1 = r, \,y_1 = \cos \theta$ gives
$$f(r) = \int_{\mathcal S} e^{-\alpha \lVert x - y \rVert^2} dS(y) =
C_n \int_0^\pi e^{-\alpha (1 + r^2 - 2 r \cos \theta)}
 \sin^{n - 2} \theta \,d\theta = \\
C_n \sqrt \pi \,\Gamma {\left( \frac {n - 1} 2 \right)}
 (\alpha r)^{-n/2 + 1} e^{-\alpha (1 + r^2)} I_{n/2 - 1}(2 \alpha r), \\
f'(r) = 2 C_n \sqrt \pi \,\Gamma {\left( \frac {n - 1} 2 \right)}
 \alpha (\alpha r)^{-n/2 + 1} e^{-\alpha (1 + r^2)}
 (I_{n/2}(2 \alpha r) - r I_{n/2 - 1}(2 \alpha r)).$$
The goal is to prove that the ratio $\omega(r) = r I_{n/2 - 1}(2 \alpha r)/I_{n/2}(2 \alpha r)$ is monotonous in $r$. Since
$$\omega'(r) = 2 \alpha r \left( 1 - \frac
 {I_{n/2 - 1}(2 \alpha r) I_{n/2 + 1}(2 \alpha r)}
 {I_{n/2}^2(2 \alpha r)} \right),$$
we need to prove
$$I_\nu^2(z) > I_{\nu - 1}(z) I_{\nu + 1}(z), \quad\nu > 0, \,z > 0.$$
Multiplying the series expansions yields
$$I_\nu^2(z) = \sum_{k \geq 0} \frac
 {\Gamma {\left( k + \nu + \frac 1 2 \right)}}
 {\sqrt \pi \,\Gamma(k + 1) \Gamma(k + \nu + 1) \Gamma(k + 2 \nu + 1)}
 z^{2 (k + \nu)}, \\
I_{\nu - 1}(z) I_{\nu + 1}(z) = \sum_{k \geq 0} \frac
 {\Gamma \!\left( k + \nu + \frac 1 2 \right)}
 {\sqrt \pi \,\Gamma(k + 1) \Gamma(k + \nu + 1) \Gamma(k + 2 \nu + 1)}
 \frac {k + \nu} {k + \nu + 1}
 z^{2 (k + \nu)}, \\
\frac
 {[z^{2 (k + \nu)}] I_\nu^2(z)}
 {[z^{2 (k + \nu)}] I_{\nu - 1}(z) I_{\nu + 1}(z)} =
1 + \frac 1 {k + \nu},$$
proving the inequality. Since $\omega(\infty) = \infty$, $\omega$ increases to infinity monotonously.
Finally, $\omega(0) = n/(2 \alpha)$. If $\omega(0) < 1$, there is a unique $r_0 > 0$ where $\omega(r_0) = 1$ and $f'(r_0) = 0$, therefore $f(r)$ increases on $[0, r_0]$ and decreases on $[r_0, \infty)$. If $\omega(0) \geq 1$, $f(r)$ decreases on $[0, \infty)$.
