Функция Грина для случайно-неоднородной среды

Главным образом интерес к вопросу распространения волн в случайно-неоднородных средах (какой является, например, атмосфера) можно объяснить бурным развитием спутниковых технологий. В этом случае становится важной задача расчета характеристик (например, амплитуды) волны, прошедшей через среду, и установления их связей с параметром неоднородности среды. Важную роль здесь играет функция Грина для случайно-неоднородной среды, зная которую можно определить эти характеристики. Рассматривается прохождение света через среду с флуктуирующей диэлектрической проницаемостью.

Волновое уравнение в случайно неоднородной среде. Функция Грина

Рассеяние электромагнитных волн в такой среде определяется системой уравнений Максвелла. Основные отличительные черты рассеяния можно рассматривать для упрощенной модели: скалярного поля $u=u(𝐫,t)$. Этим скалярным полем заменяются напряженности электрического и магнитного полей; тогда $u=u(𝐫,t)$ удовлетворяет волновому уравнению:

$\frac{\epsilon (𝐫,t)}{c^2}\frac{{\partial }^{2}u(𝐫,t)}{\partial t^2}-\mathrm{\Delta }u(𝐫,t)=0,$

$\epsilon (𝐫,t)={\epsilon }_0+\delta \epsilon (𝐫,t),$

где $c$ — скорость света в вакууме, ${\epsilon }_0$ — среднее значение диэлектрической проницаемости, $\delta \epsilon (𝐫,t)$ — флуктуации диэлектрической проницаемости. Обратим внимание, что среднее значение диэлектрической проницаемости ${\epsilon }_0$ предполагается не зависящим от координат и времени, то есть в среднем система однородна и изотропна. Также с хорошей точностью в первом приближении можно считать, что и неусреднённая диэлектрическая проницаемость $\epsilon (𝐫,t)=\epsilon (𝐫)$ не зависит от времени. Это объясняется тем, что характерные времена, отвечающие за молекулярные процессы в системе, на несколько порядков больше характерных времён электромагнитного поля, поэтому среда как бы «не успевает среагировать» на изменение поля.

Волновое уравнение с такой диэлектрической проницаемостью на самом деле является примером стохастического уравнения, так как в нём присутствует случайный параметр $\delta \epsilon (𝐫)$. Этот параметр входит в уравнение с помощью умножения, то есть мультипликативно, а не с помощью сложения (аддитивно), как в известном уравнении для броуновского движения.

Описывая рассеяние, интересны характеристики поля $u=u(𝐫,t)$, усреднённые по флуктуациям диэлектрической проницаемости. Эти характеристики: среднее значение поля $⟨u(𝐫,t)⟩$ и интенсивность $I=I(𝐫,t)$, которую определяет средний квадрат поля (усреднение так же ведётся по флуктуациям диэлектрической проницаемости) $⟨{u(𝐫,t)}^2⟩$. Статистику флуктуаций считаем заданной, а также учитываем, что усреднённое отклонение от среднего значения диэлектрической проницаемости равно нулю:

$⟨\delta \epsilon (𝐫)⟩=0.$

Начальное однородное волновое уравнение всегда имеет решение в виде $u(𝐫,t)=0$. Это очевидное тривиальное решение. Легко показать, что при отсутствии флуктуаций ненулевым решением является плоская монохроматическая волна вида:

$u(𝐫,t)=u_0\exp [i{𝐤}_{\mathrm{𝟎}}𝐫-iwt].$

Подставим это выражение в волновое уравнение. Получаем:

$-\frac{w^2{\epsilon }_0}{c^2}u(𝐫,t)+k_0^{2}u(𝐫,t)=0.$

Отсюда ясно, что предложенное решение будет удовлетворять уравнению, если частота плоской волны $w$ и волновой вектор $k_0$ связаны дисперсионным соотношением:

$k_0^2=\frac{w^2{\epsilon }_0}{c^2}.$

Понятно, что любая линейная комбинация волн, отвечающих дисперсионному соотношению, тоже является решением волнового уравнения в отсутствие флуктуаций диэлектрической проницаемости.

Определим функцию Грина $G(𝐫,t)$. Пусть эта функция является решением начального волнового уравнения, в правую часть которого добавлен расположенный в начале координат монохроматический источник (частота источника $w$). Полагаем, что источник «адиабатически включился в бесконечно далёком прошлом», для этого дополним правую часть множителем $\exp [\alpha t]$, где $\alpha$ — малая положительная величина. В окончательных выражениях будем устремлять её к нулю. Итого функция Грина удовлетворяет уравнению:

$\frac{\epsilon (𝐫)}{c^2}\frac{{\partial }^{2}G(𝐫,t)}{\partial t^2}-\mathrm{\Delta }G(𝐫,t)=e^{-iwt+\alpha t}\delta (𝐫).$

Удобно искать решение этого уравнения в виде $G(𝐫,t)=e^{-iwt+\alpha t}G(𝐫)$. Подставляя это выражение в уравнение для функции Грина, получаем:

$\frac{\epsilon (𝐫)G(𝐫)}{c^2}\frac{{\partial }^2{e}^{-iwt+\alpha t}}{\partial t^2}-e^{-iwt+\alpha t}\mathrm{\Delta }G(𝐫)=e^{-iwt+\alpha t}\delta (𝐫).$

От двойного дифференцирования экспоненты по времени появится множитель $-(w+i\alpha {)}^2$, тогда получаем уравнение для функции $G(𝐫)$:

$\frac{-\epsilon (𝐫)(w+i\alpha {)}^2}{c^2}G(𝐫)-\mathrm{\Delta }G(𝐫)=\delta (𝐫).$

Нужно решить это уравнение для некоторой диэлектрической проницаемости $\epsilon (𝐫),$ а затем усреднить это решение по всевозможным отклонениям $\delta \epsilon (𝐫)$. Но оказывается, что нет возможности получить решение этого уравнения при произвольной диэлектрической проницаемости, поэтому решение ищется с использованием теории возмущений, полагая отклонение $\delta \epsilon (𝐫)$ малой величиной.

Функция Грина для среды без флуктуаций диэлектрической проницаемости

Для начала необходимо найти функцию Грина $G(𝐫,t)$, отвечающую волновому уравнению без отклонений диэлектрической проницаемости, то есть $\epsilon (𝐫)={\epsilon }_0$: Шаблон:Формула Снова ищем решение в виде $G_0(𝐫,t)=e^{-iwt+\alpha t}{G}_0(𝐫)$. Тогда $G_0(𝐫)$ удовлетворяет уравнению: Шаблон:Формула где величиной $k_0=\frac{\sqrt{{\epsilon }_0}(w+i\alpha )}{c}$. Видно, что у $k_0$ присутствует мнимая положительная часть, далее нам это понадобится. Уравнение $(2)$ удобно решать с помощью преобразования Фурье вида: Шаблон:Формула Шаблон:Формула Выражение $(3)$ — прямое Фурье-преобразование, $F(𝐤)$ — Фурье-образ функции $F(𝐫)$, выражение $(4)$ — обратное Фурье-преобразование. Образ функции Грина $G_0(𝐫)$ будем обозначать через $G_0(𝐤)$. Применяя Фурье-преобразования к уравнению $(2)$ и учитывая, что $\delta$-функция является Фурье-образом единицы, получаем: Шаблон:Формула Шаблон:Формула Чтобы получить функцию $G_0(𝐫)$, делаем обратное Фурье-преобразование $G_0(𝐤)$: Шаблон:Формула Будем вычислять этот интеграл в сферической системе координат, выбрав полярную ось вдоль вектора $𝐫$ (под полярной осью понимается ось, от которой отсчитывается угол $\theta$):

$G_0(𝐫)=\frac{1}{(2\pi {)}^3}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{e^{i𝐤𝐫}}{k^2-k_0^2}d𝐤=\frac{1}{(2\pi {)}^3}\underset{0}{\overset{2\pi }{\int }}d\phi \underset{0}{\overset{\pi }{\int }}\sin \theta d\theta \underset{0}{\overset{\mathrm{\infty }}{\int }}{k}^2\frac{e^{ikr\cos \theta }}{k^2-k_0^2}dk=$

$=\frac{2\pi }{8{\pi }^3}\underset{\pi }{\overset{0}{\int }}d\cos \theta \underset{0}{\overset{\mathrm{\infty }}{\int }}{k}^2\frac{e^{ikr\cos \theta }}{k^2-k_0^2}dk=\frac{1}{4{\pi }^2}\underset{0}{\overset{\mathrm{\infty }}{\int }}\frac{k^2}{k^2-k_0^2}[\frac{e^{ikr\cos \theta }}{ikr\cos \theta }{]}_{\pi }^{0}dk=$

$=\frac{1}{4{\pi }^2}\underset{0}{\overset{\mathrm{\infty }}{\int }}\frac{k^2}{k^2-k_0^2}(\frac{e^{ikr}}{ikr}-\frac{e^{-ikr}}{ikr})dk=\frac{1}{4{\pi }^{2}ir}\frac{1}{2}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{k(e^{ikr}-e^{-ikr})}{k^2-k_0^2}dk=$

$=\frac{1}{8{\pi }^{2}ir}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{k{e}^{ikr}}{k^2-k_0^2}dk-\frac{1}{8{\pi }^{2}ir}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{k{e}^{-ikr}}{k^2-k_0^2}dk=\frac{2\pi i{e}^{i{k}_{0}r}}{8{\pi }^{2}ir2}+\frac{2\pi i{e}^{i{k}_{0}r}}{8{\pi }^{2}ir2}=\frac{e^{i{k}_{0}r}}{4\pi r}.$

Для вычисления интеграла по сферическим координатам мы воспользовались четностью функции $\frac{k(e^{ikr}-e^{-ikr})}{k^2-k_0^2}$, а также последние интегралы брались по вычетам. Для первого слагаемого контур интегрирования замыкался сверху, в этой полуплоскости затухает $e^{ikr}$, тогда вычет берётся в $k=k_0=\frac{\sqrt{{\epsilon }_0}(w+i\alpha )}{c}$. Для второго слагаемого замыкали контур в нижней полуплоскости, и тогда срабатывает вычет в точке $k=-k_0$, при этом необходимо не забыть, что контур обходится по часовой стрелке, тогда как в теореме по вычетам используется обход против часовой стрелки. Направление обхода можно легко изменить, добавив во втором слагаемом множитель $(-1)$.

Итоговое выражение для функции Грина будет:

$G_0(𝐫,t)=\frac{e^{i{k}_{0}r-iwt+\alpha t}}{4\pi r}.$

Это расходящаяся сферическая волна. Амплитуда этой волны убывает как $\frac{1}{r}$ по мере удаления от источника.

Функция Грина с учётом флуктуаций диэлектрической проницаемости

Перепишем уравнение

$-\frac{\epsilon (𝐫)(w+i\alpha {)}^2}{c^2}G(𝐫)-\mathrm{\Delta }G(𝐫)=\delta (𝐫)$

в виде

$(-k_0^2-\mathrm{\Delta })G(𝐫)=\frac{\delta \epsilon (𝐫)}{{\epsilon }_0}{k}_0^{2}G(𝐫)+\delta (𝐫).$

Для использования теории возмущения, в которой мы будем считать $\delta \epsilon (𝐫)$ малой величиной, удобнее перейти к интегральному аналогу предыдущего уравнения:

$G(𝐫)=G_0(𝐫)+\frac{{k_0}^2}{{\epsilon }_0}\int G_0(𝐫-{𝐫}_{\mathrm{𝟏}})\delta \epsilon ({𝐫}_{\mathrm{𝟏}})G({𝐫}_{\mathrm{𝟏}})d{r}_1.$

Тогда можно легко написать итерационное решение в виде ряда:

$G(𝐫)=G_0(𝐫)+\frac{k_0^2}{{\epsilon }_0}\int G_0(𝐫-{𝐫}_{\mathrm{𝟏}})\delta \epsilon ({𝐫}_{\mathrm{𝟏}})G_0({𝐫}_{\mathrm{𝟏}})d{𝐫}_{\mathrm{𝟏}}+$

$+\frac{k_0^4}{{\epsilon }_0^2}\int G_0(𝐫-{𝐫}_{\mathrm{𝟏}})\delta \epsilon ({𝐫}_{\mathrm{𝟏}})\int G_0({𝐫}_{\mathrm{𝟏}}-{𝐫}_{\mathrm{𝟐}})\delta \epsilon ({𝐫}_{\mathrm{𝟐}})G_0({𝐫}_{\mathrm{𝟐}})d{𝐫}_{\mathrm{𝟏}}d{𝐫}_{\mathrm{𝟐}}+\dots$

Величина $G(𝐫)$ — случайная величина. В дальнейшем её необходимо усреднить по всевозможным флуктуациям диэлектрической проницаемости. Это представляет собой следующий трудоёмкий шаг.

Литература

• Рытов С. М., Кравцов Ю. А., Татарский В. И. Введение в статистическую радиофизику. — Ч. 2. Случайные поля. — М.: Наука, 1978.
• Исимару А. Распространение и рассеяние волн в случайно-неоднородных средах. — Т. 1, 2. — М.: Мир, 1981.