I refer to this paper. 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.

I refer to this paper. 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.

I refer to this paper. 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^{-kn_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.

I refer to this paper. 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.