Лицензия MIT
© Алексей Александрович Щербаков, 2024
В данной лекции напомним известные определения, а также остановимся на фактах, которые будут использованы в курсе при построении численных решений волнового уравнения.
Для скалярного уравнения Гельмгольца \begin{equation*} (\nabla^{2}+k^2)\psi(\boldsymbol{r}) = f(\boldsymbol{r}) \end{equation*} уравнение на функцию Грина имеет вид \begin{equation*} (\nabla^{2} + k^{2})g(\boldsymbol{r} - \boldsymbol{r}') = -\delta(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} Его решение в однородной изотропной среде с условием Зоммерфельда на бесконечности в явном виде записывается как \begin{equation*} g_0(\boldsymbol{r} - \boldsymbol{r}') = \frac{\exp\left(ik\left|\boldsymbol{r} - \boldsymbol{r}'\right|\right)} {4\pi\left|\boldsymbol{r} - \boldsymbol{r}'\right|} \end{equation*} Тогда для заданного источника $f$ решение записывается в виде интеграла \begin{equation*} \psi\left(\boldsymbol{r}\right) = -\int_{V} g_0\left(\boldsymbol{r} - \boldsymbol{r}'\right) f\left(\boldsymbol{r}'\right) d^{3}\boldsymbol{r}' \end{equation*}
Для соответствующего векторного уравнения с электрическими источниками \begin{equation*} \nabla\times\nabla\times\boldsymbol{E} - k^{2}\boldsymbol{E} = i\omega\mu\boldsymbol{J}_e \end{equation*} уравнение на тензорную функцию Грина имеет вид \begin{equation*} \nabla\times\nabla\times\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') - k^{2}\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') = \mathbb{I} \delta(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} Его решение в однородной изотропной среде с условиями Сильвера-Мюллера на бесконечности выражается через скалярную функцию Грина \begin{equation*} \mathcal{G}_0(\boldsymbol{r} - \boldsymbol{r}') = \left(\mathbb{I} + \dfrac{1}{k^{2}} \nabla\nabla \right) g(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} Решение уравнения Гельмгольца во всем пространстве \begin{equation*} \boldsymbol{E} = i\omega\mu \int_{V} \mathcal{G}_0(\boldsymbol{r} - \boldsymbol{r}') \boldsymbol{J}_e \left(\boldsymbol{r}'\right) d^{3}\boldsymbol{r}' \end{equation*}
В явном виде дифференцирование скалярной функции Грина приводит к выражению \begin{equation*} \mathcal{G}_0(\boldsymbol{r} - \boldsymbol{r}') = \dfrac{e^{ikR}}{4\pi R} \left[ \mathbb{I}\left(1+\dfrac{i}{kR}-\dfrac{1}{k^{2}R^{2}}\right)+\hat{\boldsymbol{e}}_r \hat{\boldsymbol{e}}_r^T \left(\dfrac{3}{k^{2}R^{2}} - \dfrac{3i}{kR} - 1\right)\right] \end{equation*} где $R = |\boldsymbol{r}-\boldsymbol{r}'|$. Данное выражение в точке $\boldsymbol{r}$ пропорционально полю электрического диполя, расположенного в точке с координатами $\boldsymbol{r}'$.
На больших расстояниях $r\gg r'$ \begin{equation*} \mathcal{G}_{0}(\boldsymbol{r}, \boldsymbol{r}') \sim \left(\mathbb{I} - \hat{\boldsymbol{e}}_r \hat{\boldsymbol{e}}_r^T\right)\dfrac{e^{ikr}}{4\pi r}e^{-ik\hat{\boldsymbol{r}}\cdot\boldsymbol{r}'} \end{equation*} что представляет собой поперечную cферическую волну, амплитуда которой убывает обратно пропорционально расстоянию до источника. По аналогии с электрическим полем, можно написать условие на бесконечности \begin{equation*} \lim_{R\rightarrow\infty} R \left[ \nabla\times\mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}') - ik_m\hat{\boldsymbol{R}} \times \mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}') \right] = 0 \end{equation*} где $\boldsymbol{R} = \boldsymbol{r}-\boldsymbol{r}'$.
Также, из соотношения взаимности следует \begin{equation*} \mathcal{G}_{0}\left(\boldsymbol{r},\boldsymbol{r}'\right) = \mathcal{G}_{0}^{T}\left(\boldsymbol{r}',\boldsymbol{r}\right) \end{equation*}
Если разделить потенциал на однородную и неоднородную части \begin{equation*} \varepsilon = \varepsilon_m + \Delta\varepsilon\Rightarrow k^2 = k_m^2 + \Delta (k^2) \end{equation*} подставить в уравнение для функции Грина \begin{equation*} \nabla\times\nabla\times\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') - k_m^{2}\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') = \mathbb{I} \delta(\boldsymbol{r} - \boldsymbol{r}') + \Delta (k^2)\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} и воспользоваться объемным интегральным решением с источником \begin{equation*} \mathbb{J} = -\frac{i}{\omega\mu_0}\mathbb{I} \delta(\boldsymbol{r} - \boldsymbol{r}') -i\omega\Delta\varepsilon \mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} получим \begin{equation*} \mathcal{G}(\boldsymbol{r},\boldsymbol{r}') = \mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}') + \int_{V} \mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}'') U(\boldsymbol{r}'') \mathcal{G}(\boldsymbol{r}'',\boldsymbol{r}') d^{3}\boldsymbol{r}'' \end{equation*} Это уравнение называется уравнением Дайсона.
Если рассматривается конечный объём $\Omega\in\mathbb{R}^3$, интегральные решения уравения Гельмгольца можно получить из теоремы Грина. Для некоторой векторной функции $\boldsymbol{P}$ и тензорной функции $\mathcal{Q}$, заданного объема $V$ и ограничивающей его поверхности $S$ эта теорема имеет вид \begin{equation*} \int_{V} \left[ \boldsymbol{P} \cdot \nabla\times\nabla\times\mathcal{Q} - \left( \nabla\times\nabla\times\boldsymbol{P}\right) \cdot \mathcal{Q} \right] dV = -\oint_{S} \left[ \left(\hat{\boldsymbol{n}} \times\nabla\times\boldsymbol{P}\right) \cdot \mathcal{Q} + \left(\hat{\boldsymbol{n}}\times\boldsymbol{P}\right) \cdot \nabla\times\mathcal{Q} \right] dS \end{equation*} Если сделать подстановку \begin{equation*} \boldsymbol{P} = \boldsymbol{E},\thinspace \mathcal{Q} = \mathcal{G}_0(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} воспользовать уравнением Гельмгольца для электрического поля, уравнением на функцию Грина, и законом Фарадея, получим \begin{equation*} \boldsymbol{E}(\boldsymbol{r}') = ikZ \int_{V} \boldsymbol{J} \cdot \mathcal{G}_0 dV - \oint_{S} \left[ ikZ \left(\hat{\boldsymbol{n}} \times \boldsymbol{H}\right) \cdot \mathcal{G}_0 + \left(\hat{\boldsymbol{n}}\times\boldsymbol{E}\right) \cdot \nabla\times\mathcal{G}_0\right] dS \end{equation*} Здесь, как и ранее, введен импеданс свободного пространства $Z=\sqrt{\mu/\varepsilon}$.
Аналогично, для магнитного поля можно получить \begin{equation*} \boldsymbol{H}(\boldsymbol{r}') = \int_{V} \boldsymbol{J} \cdot \nabla\times\mathcal{G}_0 dV + \oint_{S} \left[i\dfrac{k_m}{Z}\left(\hat{\boldsymbol{n}}\times\boldsymbol{E}\right) \cdot \mathcal{G}_0 - \left(\hat{\boldsymbol{n}} \times \boldsymbol{H}\right)\cdot\nabla\times\mathcal{G}_0 \right] dS \end{equation*}
Из явного выражения для тензорной функции Грина свободного пространства следует, что в точке $\boldsymbol{r} = \boldsymbol{r}'$ функция имеет полюс третьего порядка, и объемный интеграл для поля в области источников вычисляется в смысле главного значения. Чтобы вычислить этот интеграл выделим интеграл по бесконечно малому объему вблизи полюса и воспользуемся законом сохранения заряда: \begin{split} &\nabla\nabla\int_{V} g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' = \\ &= \lim_{V_{\delta}\rightarrow0} \left[ \nabla\nabla\int_{V \backslash V_{\delta}} g_0 (\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' + \nabla\nabla\int_{V_{\delta}} g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' \right] = \\ &= \lim_{V_{\delta}\rightarrow0} \left[ \nabla\nabla\int_{V \backslash V_{\delta}} g_0 (\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' - \nabla\int_{V_{\delta}} \nabla' g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' \right] = \\ &= \lim_{V_{\delta}\rightarrow0} \left[ \nabla\nabla\int_{V \backslash V_{\delta}} g_0 (\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' - \nabla \left( \int_{S_{\delta}} g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{n} \cdot \boldsymbol{J} (\boldsymbol{r}') dS' - i\omega \int_{V_{\delta}} g_{0}(\boldsymbol{r},\boldsymbol{r}') \rho_e(\boldsymbol{r}') d^{3}\boldsymbol{r}' \right) \right] \end{split} Предполагая, что плотность заряда $\rho_e(\boldsymbol{r}')$ является непрерывной функцией координат, получаем, что последнее слагаемое равно нулю. Произведение $\boldsymbol{n} \cdot \boldsymbol{J}$ задает поверхностный заряд на поверхности $S_{\delta}$, поэтому предпоследнее слагаемое сводится к произведению постоянного множителя, зависящего от формы этой поверхности на ток точке $\boldsymbol{r}$, так что в итоге \begin{equation*} \nabla\nabla\int_{V} g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' = \lim_{V_{\delta}\rightarrow0} \nabla\nabla\int_{V \backslash V_{\delta}} g_0 (\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' - \mathcal{L} \boldsymbol{J} (\boldsymbol{r}) \end{equation*} Полностью уравнение для поля запишется как \begin{equation*} \boldsymbol{E} = i\omega\mu_{0} P.V.\int_{V} \mathcal{G}_{0} (\boldsymbol{r} - \boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3}\boldsymbol{r}' + \dfrac{1}{i\omega\varepsilon_{m}} \mathcal{L} \cdot \boldsymbol{J} (\boldsymbol{r}) \end{equation*} Например, для сферического объема \begin{equation*} \mathcal{L}_{sph} = \frac{1}{3} \mathbb{I} \end{equation*}
Рассмотрим уравнение для тензорной функции Грина. Трехмерное преобразование Фурье по переменной $\boldsymbol{r}$ даёт \begin{equation*} -\boldsymbol{k} \times \boldsymbol{k} \times \mathcal{G}_0 (\boldsymbol{k},\boldsymbol{r}') - k_{0}^{2} \mathcal{G}_0 (\boldsymbol{k},\boldsymbol{r}') = \mathbb{I} e^{-i\boldsymbol{k}\boldsymbol{r}'} \Rightarrow \mathcal{G}_0 (\boldsymbol{k},\boldsymbol{r}') = \dfrac{\mathbb{I} k_{0}^{2} - \boldsymbol{k}\boldsymbol{k}^T}{k_{0}^{2}\left(k^{2}-k_{0}^{2}\right)} e^{-i\boldsymbol{k}\boldsymbol{r}'} \end{equation*} С помощью обратного преобразования Фурье получаем спектральное разложение \begin{equation*} \mathcal{G}_0 (\boldsymbol{r},\boldsymbol{r}') = \dfrac{1}{\left(2\pi\right)^{3}} \int_{-\infty}^{\infty} d^{3}\boldsymbol{k} e^{i\boldsymbol{k}(\boldsymbol{r}-\boldsymbol{r}')} \dfrac{\mathbb{I}k_{0}^{2} - \boldsymbol{k}\boldsymbol{k}^T} {k_{0}^{2}(k^{2}-k_{0}^{2})} \end{equation*} Этот интеграл расходится в классическом смысле - дробь под интегралом стремится к постоянной при $|\boldsymbol{k}|\rightarrow\infty$. Рассмотрим предел $|k_z|\rightarrow\infty$ при ограниченных $|k_{x,y}|<const$: \begin{equation*} \dfrac{\mathbb{I}k_{0}^{2} - \boldsymbol{k}\boldsymbol{k}^T} {k_{0}^{2}(k^{2}-k_{0}^{2})} \sim -\hat{\boldsymbol{e}}_z \hat{\boldsymbol{e}}_z^T \end{equation*} Выделим в интеграле сингулярную часть, а регулярную проинтегрируем по $k_z$ с помощью леммы Жордана, вычисляя вычеты в полюсах $\pm\sqrt{k_{0}^{2}-k_{x}^{2}-k_{y}^{2}}$: \begin{align*} \mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}') =& \dfrac{1}{\left(2\pi\right)^{3}} \int_{-\infty}^{\infty} d^{3}\boldsymbol{k} e^{i\boldsymbol{k}(\boldsymbol{r}-\boldsymbol{r}')} \left[ \dfrac{\mathbb{I}k_{0}^{2} - \boldsymbol{k}\boldsymbol{k}} {k_{0}^{2}\left(k^{2}-k_{0}^{2}\right)} + \dfrac{\hat{\boldsymbol{e}}_z \hat{\boldsymbol{e}}_z}{k_{0}^{2}} \right] - \dfrac{\hat{\boldsymbol{e}}_z \hat{\boldsymbol{e}}_z}{\left(2\pi\right)^{3}k_{0}^{2}} \int_{-\infty}^{\infty} d^{3}\boldsymbol{k} e^{i\boldsymbol{k}\left(\boldsymbol{r}-\boldsymbol{r}'\right)} = \\ =& \dfrac{i}{8\pi^{2}} \int_{-\infty}^{\infty} dk_xdk_y e^{i\boldsymbol{\varkappa} (\boldsymbol{\rho} - \boldsymbol{\rho}') + ik_{z} \left|z-z'\right|} \dfrac{\mathbb{I}k_{0}^{2} - \boldsymbol{k}_{0}\boldsymbol{k}_{0}}{k_{0}^{2}k_{z}} - \dfrac{\hat{\boldsymbol{e}}_z \hat{\boldsymbol{e}}_z}{k_{0}^{2}}\delta\left(\boldsymbol{r}-\boldsymbol{r}'\right) \end{align*} где $\boldsymbol{\varkappa} = (k_x,k_y)^T$, $\boldsymbol{\rho} = (x,y)^T$. Полученное выражение есть разложение тензорной функции Грина свободного пространства по плоским волнам. Подынтегральное выражение содержит только поперечные плоские волны, а сингулярность возникает только в области источника. Комопненты волнового вектора удовлетворяют дисперсионному уравнению, а знак третьей компоненты выбирается исходя из условий излучения на бесконечности: \begin{equation*} \boldsymbol{k}_{0} = \begin{cases} \hat{\boldsymbol{e}}_x k_{x} + \hat{\boldsymbol{e}}_y k_{y} + \hat{\boldsymbol{e}}_z k_{z} & z-z'>0 \\ \hat{\boldsymbol{e}}_x k_{x} + \hat{\boldsymbol{e}}_y k_{y} - \hat{\boldsymbol{e}}_z k_{z} & z-z'<0 \end{cases},\thinspace\Im m\left[k_{z}\right] \geq 0 \end{equation*}
Ранее мы получили полный базис векторных сферических волн $\boldsymbol{M}$, $\boldsymbol{N}$, $\boldsymbol{L}$, первые две из которых являются соленодальными, а третья - безвихревой. В декартовых координатах произвольное поле может быть разложено по этим функциям как \begin{equation*} \boldsymbol{E} = \iiint_{-\infty}^{\infty} d^{3}\boldsymbol{k} \left[ a_M(\boldsymbol{k})\boldsymbol{M}(\boldsymbol{k},\boldsymbol{r}) + a_N(\boldsymbol{k})\boldsymbol{N}(\boldsymbol{k},\boldsymbol{r}) + a_L(\boldsymbol{k})\boldsymbol{L}(\boldsymbol{k},\boldsymbol{r}) \right] \end{equation*} Учитывая соотношения ортогональности, можно записать разложение единицы (его можно проверить прямой подстановкой) \begin{equation*} \mathcal{I}\delta\left(\boldsymbol{r} - \boldsymbol{r}'\right) = \dfrac{1}{\left(2\pi\right)^{3}} \iiint_{-\infty}^{\infty} d^{3}\boldsymbol{k} \left[ \dfrac{\boldsymbol{M}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{M}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right) + \boldsymbol{N}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{N}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)}{\varkappa^{2}} + \dfrac{\boldsymbol{L}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{L}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)} {k_0^{2}}\right] \end{equation*} Чтобы выразить функцию Грина, запишем для нее разложение, аналогичное разложению поля по векторным волнам: \begin{equation*} \mathcal{G}=\iiint_{-\infty}^{\infty} d^{3}\boldsymbol{k} \left[ \boldsymbol{M}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{a}_M^T\left(\boldsymbol{k},\boldsymbol{r}'\right) + \boldsymbol{N}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{a}_N^T\left(\boldsymbol{k},\boldsymbol{r}'\right) + \boldsymbol{L}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{a}_L^T\left(\boldsymbol{k},\boldsymbol{r}'\right) \right] \end{equation*} Подстановка этого разложения и разложения единицы в уравнение Гельмгольца для функции Грина, и использование соотношений ортогональности дает явное выражение коэффициентов, так что для тензора Грина имеем \begin{equation*} \mathcal{G} = \dfrac{1}{\left(2\pi\right)^{3}} \iiint_{-\infty}^{\infty} d^{3}\boldsymbol{k} \left[ \dfrac{\boldsymbol{M}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{M}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)+\boldsymbol{N}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{N}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)}{\left(k^{2}-k_{0}^{2}\right)\varkappa^{2}} - \dfrac{\boldsymbol{L}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{L}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)}{k_{0}^{2}k^{2}}\right] \end{equation*} По аналогии со спектральным разложением, выделение сингулярной части и интегрирование с применением леммы Жордана приводит к разложению тензора Грина по плоским волнам, эквивалентного полученному выше при рассмотрении спектрального разложения: \begin{equation*} \mathcal{G} = \dfrac{i}{8\pi^{2}}\iint_{-\infty}^{\infty} dk_xdk_y \left[ \dfrac{\boldsymbol{M}\left(\boldsymbol{\varkappa},\pm k_{z},\boldsymbol{r}\right)\boldsymbol{M}^T\left(-\boldsymbol{\varkappa},\mp k_{z},\boldsymbol{r}\right) + \boldsymbol{N}\left(\boldsymbol{\varkappa},\pm k_{z},\boldsymbol{r}\right)\boldsymbol{N}^T\left(-\boldsymbol{\varkappa},\mp k_{z},\boldsymbol{r}\right)}{k_{z}\varkappa^{2}}\right] -\dfrac{\hat{\boldsymbol{e}}_z\hat{\boldsymbol{e}}_z^T}{k_{0}^{2}}\delta\left(\boldsymbol{r}-\boldsymbol{r}'\right) \end{equation*} Здесь также подразумевается, что компоненты волнового вектора удовлетворяют дисперсионному уравнению и применяется условие на выбор знака перед константой распространения вдоль оси $Z$.
Аналогично случаю плоских волн получается выражение для тензорной функции Грина, разложенной по векторным сферическим волнам: \begin{equation*} \mathcal{G}\left(\boldsymbol{r},\boldsymbol{r}'\right) = ik_{0}\sum_{nm}\dfrac{1}{n\left(n+1\right)}\left[\boldsymbol{M}^{\sigma}_{nm}\left(k_{m}\boldsymbol{r}\right)(\hat{\boldsymbol{M}}^{\sigma'}_{n,-m})^T\left(k_{m}\boldsymbol{r}'\right)+\boldsymbol{N}^{\sigma}_{nm}\left(k_{m}\boldsymbol{r}\right)(\hat{\boldsymbol{N}}^{\sigma'}_{n,-m})^T\left(k_{m}\boldsymbol{r}'\right)\right]-\dfrac{\hat{\boldsymbol{e}}_r\hat{\boldsymbol{e}}_r}{k_{0}^{2}}\delta\left(\boldsymbol{r}-\boldsymbol{r}'\right) \end{equation*} Здесь для двух типов решений $\sigma=1$, $\sigma'=3$ при $r<r'$, и $\sigma=3$, $\sigma'=1$ при $r>r'$.