Given a spherical reference frame $\left(\rho,\phi,\theta\right)$, the reduced wave equation can be written as:
$$\nabla^2U=k^2n^2U$$
in which:
$U=U(\rho,\phi,\theta)$
The solutions of this equation can be found using classical methods of PDE knowing the initial and boundary conditions.
Suppose we have, $k$ constant with $k\in\mathbb R$:
$n=n\_0+\left(\nu\_{\rho}+\nu\_{\phi}+\nu\_{\theta}\right)$ where 
$\nu\_{\rho},\nu\_{\phi},\nu\_{\theta}$ are normally distribuited functions with zero mean and $\sigma(\nu\_{\rho})=\sigma\_{\rho}, \sigma(\nu\_{\phi})=\sigma\_{\phi},\sigma(\nu\_{\theta})=\sigma\_{\theta}$, how can be found the statistical properties of $U(\rho,\phi,\theta)$? Thanks.



Solutions of a stochastic reduced wave equation
I refer to [this paper](https://www.researchgate.net/publication/228448963_Stochastic_Partial_Differential_Equations_as_priors_in_ensemble_methods_for_solving_inverse_problems). There is a straightforward approach in this case. You can separate the refraction index in the following way
$$
 n^2=n\_0^2+\xi
$$
where $\xi$ is the random part. Then, the equation is
$$
 \nabla^2U+k^2n\_0^2U=-\xi k^2U.
$$
Now, turning to the paper I just cited, take the Green function to be ($d=3$ for the sake of simplicity)
$$
 G(|{\bf r}|) =\frac{e^{ikn\_0|{\bf r}|}}{4\pi |{\bf r}|}
$$
and the differential equation turns into an integral one
$$
 U({\bf r}) = U\_0({\bf r}) - k^2\int d^3x'G(|{\bf r}-{\bf r'}|)\xi({\bf r'})U({\bf r'}).
$$
One can iterate this equation and obtain the following solution series
$$
 U({\bf r}) = U\_0({\bf r}) - k^2\int d^3x'G(|{\bf r}-{\bf r'}|)\xi({\bf r'})U\_0({\bf r'})+k^4\int d^3x'G(|{\bf r}-{\bf r'}|)\xi({\bf r'})U\_0({\bf r'})\times
\int d^3x''G(|{\bf r'}-{\bf r''}|)\xi({\bf r''})U\_0({\bf r''})+\ldots.
$$
By averaging on the distribution for $\xi$, a normal one by hypothesis, one gets a solution series in terms of the correlation function of the random variable. A simple case is this
$$
 \langle\xi({\bf r})\rangle =0 \qquad \langle\xi({\bf r})\xi({\bf r'})\rangle=\xi\_0^2\delta^3({\bf r}-{\bf r'})
$$
yielding
$$
 U({\bf r}) \stackrel{?}{=} U\_0({\bf r})+k^4\xi\_0^2G(0)\int d^3x'G(|{\bf r}-{\bf r'}|)U\_0^2({\bf r'})+\ldots.
$$
that is not well defined mathematically. You can consider a case with a finite volume and recover from this situation.


