I want to evaluate the following surface integral
$$\int_S d S \, \, \exp\left(\frac{\Gamma^2 (x^2 \sigma_x^2 + y^2 \sigma_y^2 + z^2 \sigma_z^2) - 2 c R \Gamma (x x_0 + y y_0 + z z_0)}{2 c^2 R^2}\right)$$
where $S$ describes the surface of a sphere of radius $R$, the real positive parameters $\Gamma, x_0, y_0, z_0$ define the physical system being considered, $c$ is the speed of light. The general point $\bf{r}$ = $(x, y, z)$ lies on the surface, $S$, of the sphere with radius $R$ that I am trying to integrate over.
How can I obtain an analytic expression for this? My current attempts have made no progress. I have tried Mathematica. I have also changed the integral into spherical coordinates using $d S = R^2 \sin\theta\ d\theta\ d\phi$ and with $x = R\sin\theta\cos\phi, y = R\sin\theta\sin\phi, z = R\cos\theta$. However, I am having difficulties in evaluating the integrals which result from this change in basis as well. For example if we consider the integral wrt $\phi$ from 0 to $2 \pi$ which results, we have:
$$\int_0^{2\pi}d\phi \, \, \exp\left(\frac{\Gamma^2 (\cos^2\phi\sin^2\theta \sigma_x^2 + \sin^2\phi\sin^2\theta \sigma_y^2) - 2 c \Gamma (x_0\cos\phi\sin\theta + y_0\sin\phi\sin\theta)}{2 c^2}\right)$$
This appears to no simpler an evaluation.
Are there any other methods or substitutions that can be used to yield an analytical expression for the first integral written above?
ADDENDUM: PHYSICAL CONTEXT OF PROBLEM: As a final note, the physical system that I am considering is the propagation of a photon to all points along the surface of a sphere. That is along some source plane, I have a single photon that has a 3D-Gaussian distribution (or [Multivariate][2] with 0 off diagonal terms in the [covariance matrix][3]).
As a result I have a ket (state of the photon), which is proportional to this 3D-Gaussian. In fact you can 'see' the resemblence of a Gaussian as the exponential term in the above integral. The difference is due to the propagation of this source photon to the surface of a sphere a distance $R$ away from the source. The propagator that I use is the retarded Green's function from electrostatics.
The requirement of normalising the state along the surface of a sphere (for unit probablities) encounters the integral above.
The following assumes that the $\sigma_x$, $\sigma_y$, and $\sigma_z$ are equal to a common value $\sigma$ - doing so allows one to make use of spherical symmetry.
Write $S_R$ for the sphere of radius $R$. Then, the above integral can be labeled and rewritten as $$ \phi_{(R, \Gamma, \sigma )}(y) = \int_{S_R} d S_R \, \, \exp\left(\frac{(\Gamma\sigma)^2 R^2 - 2 c R \Gamma \,x\cdot y}{2 c^2 R^2}\right),$$ where
Taking out the constant factor and simplifying, one gets $$ \phi_{(R, \Gamma, \sigma )}(y) = \exp \left(\,(\Gamma\sigma)^2/2 c^2 \right)\int_{S_R} d S_R \, \, \exp\left(\frac{ - \Gamma \,x\cdot y}{ c R}\right).$$
For neatness, group the terms $y$, $\Gamma$ and $c$ of the integral together, and replace the integral with one over the unit sphere (so $dS_R = R^2 dS_1$, and the $R$ cancels in the exponent, as we are replacing $x$ with $Rx$):
$$ \phi_{(R, \Gamma, \sigma )}(y) = \exp \left(\,(\Gamma\sigma)^2/2 c^2 \right)\, R^2 \,\phi (\Gamma y/c),$$ where
$$ \phi( y ) = \int_{S_1} d S_1 \, \, \exp\,\left( - x\cdot y\right).$$
Now, by symmetry, $\phi (y)$ only depends on the magnitude of $y$, so $\phi (y) = f( \|y\| )$, where $$ f(r ) = \int_{S_1} d S_1 \, \, \exp\, \left(- x\cdot kr\right),$$ with $r\ge 0$ and $k=(0,0,1)$.
Now, $f(0) = 4 \pi$.
For $r\not = 0$:
If $\theta$ is the angle between the $k$ and $x$ ($0$ pointing northwards, $\pi$ pointing southwards) we can rewrite the above as $$f(r) = 2\pi\int_0^\pi \exp\, \left(- r\cos \theta \right) \,\sin\theta \,d\theta.$$ Finally $$ f(r) = 2 \pi \, \left( e^r - e^ {-r} \over r \right) = {4 \pi \sinh r \over r}.$$