Reflection from the periodically modulated interface of a hyperbolic metamaterial

TL;DR

This paper numerically investigates how periodically modulated hyperbolic metamaterial interfaces can minimize reflection and induce subwavelength field localization, providing analytical and numerical insights into their scattering properties.

Contribution

It introduces an optimal sawtooth-shaped interface profile for minimal reflection and generalizes a numerical method using Green's second identity for field calculations.

Findings

Sawtooth interface profile minimizes reflection

Subwavelength localization occurs at specific orientations

Analytical eigenmode expression derived

Abstract

A numerical study of the features of the scattering of incident radiation at a periodically modulated hyperbolic metamaterial boundary was carried out. We have found optimal - sawtooth-shaped profile of the interface which provides minimal reflection. It is shown that at a certain orientation of the optical axes of the metamaterial, a subwave localization of the field occurs at the saw-tooth interface. An approximate analytical expression for the eigenmode was found. A generalization of the numerical method is implemented for calculating fields, based on the application of the Green's second identity.

Equations

$$x_B=\pi \int \cos \theta d\tau -P\sin (2\tau ),y_B=\pi \int \sin \theta d\tau ,$$

$${\epsilon }_{\Vert }={\epsilon }_M{\rho }_M+{\epsilon }_D(1-{\rho }_M),$$

$${\epsilon }_{\perp }={\epsilon }_D\frac{{\epsilon }_M(1+{\rho }_M)+{\epsilon }_D(1-{\rho }_M)}{{\epsilon }_D(1+{\rho }_M)+{\epsilon }_M(1-{\rho }_M)},$$

$${\epsilon }_{\varphi \varphi }^{-1}\frac{1}{r}\frac{\partial }{\partial r}r\frac{\partial h}{\partial r}+{\epsilon }_{rr}^{-1}\frac{1}{r^2}\frac{{\partial }^{2}h}{\partial {\varphi }^2}+k_0^{2}h=0.$$

$$\frac{d^{2}H}{d{\zeta }^2}-\frac{{\epsilon }_{\varphi \varphi }}{{\epsilon }_{rr}}{\left(\frac{\pi }{2{\varphi }_0}\right)}^{2}H+k_0^2{\epsilon }_{\varphi \varphi }{e}^{2\zeta }H=0.$$

$$k_0{r}_{max}\approx {\epsilon }_{rr}^{-1/2}\frac{\pi }{2{\varphi }_0},$$

$$H\approx \sin \left(\sqrt{-{\epsilon }_{\varphi \varphi }/{\epsilon }_{rr}}\frac{\pi }{2{\varphi }_0}\log r\right).$$

$${\int }_U\left(\psi \Delta \varphi -\varphi \Delta \psi \right)dV={\oint }_{\partial U}\left(\psi \frac{\partial \varphi }{\partial \mathbf{n}}-\varphi \frac{\partial \psi }{\partial \mathbf{n}}\right)dS.$$

$$\Delta h+k_0^{2}h=0,$$

$$\Delta G+k_0^{2}G=\delta (x-x_0)\delta (y-y_0)\equiv \delta (\mathbf{r}-{\mathbf{r}}_0),$$

$$\frac{\partial G_{ij}}{\partial n_j}{h}_j^{sct}d{\tau }_j-G_{ij}\frac{\partial h_j^{sct}}{\partial n_j}d{\tau }_j=0,$$

$$\frac{1}{{\epsilon }_{YY}}\frac{{\partial }^2{h}^{tr}}{\partial X^2}+\frac{1}{{\epsilon }_{XX}}\frac{{\partial }^2{h}^{tr}}{\partial Y^2}+k_0^2{h}^{tr}=0,$$

$$\xi =\sqrt{{\epsilon }_{YY}}X,\eta =\sqrt{{\epsilon }_{XX}}Y$$

$$\frac{\partial {\tilde{G}}_{ij}}{\partial {\tilde{n}}_j}d{\tilde{\tau }}_j=\frac{{\tilde{G}}_{ij}^{\prime }{k}_0}{{\tilde{\rho }}_{ij}}[({\xi }_j-{\xi }_i)d{\eta }_j-({\eta }_j-{\eta }_i)d{\xi }_j]=$$

$$\sqrt{{\epsilon }_{XX}{\epsilon }_{YY}}\frac{{\tilde{G}}_{ij}^{\prime }{k}_0}{{\tilde{\rho }}_{ij}}[(x_j-x_i)d{y}_j-(y_j-y_i)d{x}_j]=$$

$$\sqrt{{\epsilon }_{\Vert }{\epsilon }_{\perp }}\frac{{\tilde{G}}_{ij}^{\prime }{k}_0}{{\tilde{\rho }}_{ij}}({\mathbf{r}}_j-{\mathbf{r}}_i){\mathbf{n}}_{j}d{\tau }_j,$$

$$\frac{\partial h_j^{tr}}{\partial {\tilde{n}}_j}d{\tilde{\tau }}_j=\frac{\partial h_j^{tr}}{\partial \xi }d{\eta }_j-\frac{\partial h_j^{tr}}{\partial \eta }d{\xi }_j=$$

$$\sqrt{{\epsilon }_{\Vert }{\epsilon }_{\perp }}\left(\frac{\partial h_j^{tr}}{\partial X}d{Y}_j\frac{1}{{\epsilon }_{YY}}-\frac{\partial h_j^{tr}}{\partial Y}d{X}_j\frac{1}{{\epsilon }_{XX}}\right)=$$

$$\sqrt{{\epsilon }_{XX}{\epsilon }_{YY}}i{k}_0\left(E_X^{tr}d{X}_j+E_Y^{tr}d{Y}_j\right)=\sqrt{{\epsilon }_{\Vert }{\epsilon }_{\perp }}i{k}_0{E}_{\tau }^{tr}d{\tau }_j,$$

$$\sum _j\left[\frac{{\tilde{G}}_{ij}^{\prime }{k}_0}{{\tilde{\rho }}_{ij}}({\mathbf{r}}_j-{\mathbf{r}}_i){\mathbf{n}}_j{h}_j^{tr}d{\tau }_j-i{k}_0{\tilde{G}}_{ij}(E_{\tau }^{tr}{)}_{j}d{\tau }_j\right]=0,$$

$$\sum _j\left[\frac{G_{ij}^{\prime }{k}_0}{{\rho }_{ij}}({\mathbf{r}}_j-{\mathbf{r}}_i){\mathbf{n}}_j{h}_j^{sct}d{\tau }_j-i{k}_0{G}_{ij}(E_{\tau }^{sct}{)}_{j}d{\tau }_j\right]=0.$$

$$h_j^{tr}=h_j^{sct}+h_j^{inc};(E_{\tau }^{tr}{)}_j=(E_{\tau }^{sct}{)}_j+(E_{\tau }^{inc}{)}_j.$$

$$V_i{G}_{ij}=h_j^{sct},$$

$${\tilde{V}}_i{\tilde{G}}_{ij}=h_j^{tr},$$

$$h(x,y)=V_{i}G(k_0\sqrt{(x-x_i{)}^2+(y-y_i{)}^2})+h^{inc}(x,y),$$

$$\tilde{h}(X,Y)={\tilde{V}}_{i}G(k_0\sqrt{{\epsilon }_{YY}(X-X_i{)}^2+{\epsilon }_{XX}(Y-Y_i{)}^2}),$$

Full text

I Введение

Большой интерес исследователей к так называемым гиперболическим метаматериалам (ГМ) h1 ; h2 ; h3 ; h4 ; Hashemi_PRB_12 ; Nefedov_SciRep_13 ; Nefedov_JOpt_13 обусловлен их необычными и полезными для приложений электромагнитными свойствами. ГМ это анизотропная резонансная среда, которая характеризуется изочастотной поверхностью гиперболического типа, а не эллипсоидальной, как для обычных сред. В оптике ГМ может быть реализован, например, как металло-диэлектрическая планарная наноструктура, либо как решетка металлических нанопроволок в диэлектрической матрице.

Уникальные свойства ГМ могут найти множество приложений, среди которых отрицательная рефракция 23 ; 24 , создание субволновых изображений Jacob_OE_06 ; Sun_NatCom_15 ; Tumkur_APL_12 ; Tumkur1_APL_12 ; 25 ; Mantel_OL_11 , локализация света в широком частотном диапазоне cui_NanoLett12 ; 167 . Очень высокая плотность фотонных состояний в ГМ Jacob_ApplPhysB_10 ; Shekhar_NanoConvergence_14 ; we4 обеспечивает эффективный излучательный теплообмен h4 ; Hashemi_PRB_12 ; Nefedov_SciRep_13 , что может использоваться в дизайне наноразмерных устройств нагрева/охлаждения.

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

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

II Постановка задачи

В простейшей постановке задача сводится к отражению плоской волны TM поляризации от двумерного профиля $\hat{\varepsilon}$. Метаматериал описывается диэлектрическим тензором $\hat{\varepsilon}$, и в главных осях $\hat{\varepsilon}=[\varepsilon_{\|},\ 0;\ 0,\ \varepsilon_{\bot}]$. Граница $x_{B}(\tau),y_{B}(\tau)$, разделяющая свободное пространство и метаматериал, является периодической функцией параметра $\tau$ (длина вдоль границы) и характеризуется масштабом вдоль границы $L_{x}$, амплитудой профиля $L_{y}$ и формой. Масштабирование профиля по $x$ и $y$ координатам для заданной формы выполняется элементарно, но задача выбора формы не так тривиальна. Оптимальным представляется задание $x_{B}(\tau),y_{B}(\tau)$ математической формулой с параметром, изменение которого существенно меняет форму границы. В дальнейшем мы использовали функцию, задающую $x_{B}(\tau),y_{B}(\tau)$ на длине элементарной ячейки $\tau=[-\pi/2:3/2\pi]$, в виде

[TABLE]

где $\theta=-0.4\pi\cos\tau$. Управляющий параметр $P$ меняет форму профиля от почти пилообразной к почти косинусоидальной и до напоминающей меандр. На рисунке 1 приведены в качестве иллюстрации кривые $y_{B}(x_{B})$ для набора значений $P$.

Перебирая различные значения параметра $P$, масштаба элементарной ячейки вдоль границы $L_{x}$ и амплитуды профиля $L_{y}$, можно найти оптимальное сочетание, которое обеспечивает слабое отражение падающего излучения от метаматериала. Дополнительным параметром оптимизации может служить также угол ориентации главных (оптических) осей метаматериала по отношению к осям $x$, $y$.

Коэффициент отражения от полупространства ГМ с профилированной границей может быть найден при решении уравнения Гельмгольца методом конечных элементов (FEM). Однако попытка использовать стандартную FEM программу решения уравнений в частных производных оказалась неудачной. Поскольку в этом случае нет прямого доступа к коду, то возникают проблемы: появляется отражение излучения от границ расчетного интервала, приходится рассматривать отражение не плоской волны, но широкого гауссова пучка из-за того, что во встроенном модуле нет опции задания периодических по $x$ граничных условий. Это все вносит ошибки в расчет коэффициента отражения. Более удобным оказался поэтому другой метод расчета, основанный на использовании второго тождества Грина (Green’s second identity), подробное описание которого дается в Приложении.

III Результаты численного исследования

Тестирование метода и соответствующей расчетной программы проводилось для предельных случаев (i) слабой модуляции границы, (ii) одинаковых диэлектрических характеристик двух сред, (iii) малого характерного периода модуляции, когда работает приближение эффективной среды. Расчетные и теоретические значения коэффициента отражения совпали в этих случаях с точностью $\sim 1$% (ошибка обусловлена дискретизацией вычислений). Кроме того, мы проводили сравнение распределения полей, полученного нашим методом и рассчитанного стандартной программой, которая использует метод конечных элементов. С учетом сказанного выше относительно трудностей применения стандартной программы, это сравнение показало разумное совпадение.

Параметры гиперболической среды, использованные в численном моделировании были найдены для метаматериала представляющего собой диэлектрическую матрицу с металлическими включениями в форме проволок (wire medium). В приближении эффективной среды компоненты диэлектрического тензора $\hat{\varepsilon}$ вдоль проволок $\varepsilon_{\|}$ и поперек них $\varepsilon_{\bot}$ находятся как wire_medium

[TABLE]

где $\varepsilon_{M}$ диэлектрическая проницаемость металла, $\varepsilon_{D}$ проницаемость диэлектрической матрицы и $\rho_{M}$ объемная доля металла в составе композита. Далее мы используем для расчетов параметры стекла ($\varepsilon_{D}$=2.1590) и серебра с объемной долей $\rho_{M}=0.25$. Дисперсионная зависимость проницаемости серебра $\varepsilon_{M}(\lambda)$ берется из базы данных database . Для вакуумной длины волны излучения $\lambda=0.6\mu$m величина $\varepsilon_{M}=-15.9822+0.5899i$, и расчеты по формулам (III) дают значения $\varepsilon_{\|}=4.2660+0.0318i$ и $\varepsilon_{\bot}=-2.3763+0.1475i$.

Результаты расчета рассеяния излучения при нормальном падении на периодическую границу ГМ с помощью метода, описанного в Приложении, иллюстрируются на рисунке 2, где приведены значения коэффициента отражения как функции амплитуды модуляции $A/\lambda$ и нормированного волнового числа модуляции $\varkappa/k_{0}=L_{x}/\lambda$ для двух различных типов профиля, пилообразного (с параметром формы $P=-1.2$) и косинусоидального (с параметром формы $P=0.5$).

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

На рисунке 3 приведены рассчитанные распределения поля для трех различных типов профиля границы гиперболической среды. Коэффициент отражения плоской волны $h^{inc}=\exp(-ik_{0}y)$ оказывается равным $R=$0.016 (панель (a)), 0.065 (панель (b)) и 0.16 (панель (c)). Оптимальной формой границы, обеспечивающей минимум отражения, является пилообразный профиль (см. также рис. 2). В случае (b) (косинусоидальный профиль) отражение примерно вдвое меньше, чем от плоской границы полупространства с $\varepsilon=\varepsilon_{\|}$, но максимальное отражение наблюдается от меандроподобного профиля (панель(с)), и коэффициент отражения в этом случае заметно больше, чем от плоской границы.

Коэффициент отражения существенно зависит от ориентации осей анизотропии метаматериала по отношению к границе. При повороте осей на 90o, когда компонента тензора проницаемости вдоль границы отрицательна, $\varepsilon_{xx}=\varepsilon_{\bot}=-2.3763+0.1475i$, коэффициент отражения от плоской границы близок к единице, и профилирование границы уменьшает отражение. На рисунке 4 приведена зависимость коэффициента отражения от профилированной границы с параметром формы $P=-1.2$(a) и $P=0.5$ (b) как функция амплитуды модуляции $A/\lambda$ и нормированного волнового числа модуляции $\varkappa/k_{0}=L_{x}/\lambda$. Существенное уменьшение отражения наблюдается при условии роста как амплитуды модуляции $A$, так и волнового числа модуляции $\varkappa$ для $P=-1.2$. Эта же тенденция, хотя и с гораздо меньшим эффектом, прослеживается для профиля другой формы с параметром $P=0.5$ (форма профиля изображена на рис. 3 (b)).

Следует отметить, что основным параметром, определяющим отражение (по крайней мере, при достаточно сильной модуляции границы) является отношение поперечного (амплитуды модуляции $A$) и продольного ($L_{x}=2\pi/\varkappa$) масштабов модуляции, то есть параметр $g=A\varkappa/2\pi$. На рисунке 5 приведена зависимость $R(g)$, где красными кружками (черными звездочками) отмечены данные для пилообразного, $P=-1.2$ (косинусоидального, $P=0.5$) профиля, и малая степень дисперсии коэффициента отражения при $g>3.5$ указывает на то, что одинаковый эффект достигается как при сильной (относительно длинноволновой), так и слабой (но мелкомасштабной) модуляции.

IV Субволновая собственная мода

Интересный эффект наблюдается для такой ориентации оптических осей ГМ при рассеянии падающего излучения на пилообразом профиле. При достаточно узких зубцах ‘‘пилы’’ поле локализуется вблизи ее вершины. Характерные распределения поля приведены на рисунке 6.

Эти локализованные распределения имеют субволновой характерный масштаб, и анализ показывает, что теоретически масштаб локализации может быть сколь угодно малым.

Действительно, в этом случае задача имеет приближенное аналитическое решение. Для анализа перейдем в полярную систему координат ($r$, $\phi$), с центром в вершине клина (зуба пилы) и будем считать, что диэлектрический тензор также имеет полярную симметрию, так что диагональные компоненты $\varepsilon_{\phi\phi}=\varepsilon_{\bot}<0$ и $\varepsilon_{rr}=\varepsilon_{\|}>0$. Уравнение Гельмгольца для магнитного поля в этом случае имеет вид

[TABLE]

Очевидно, что субволновая локализация поля предполагает его экспоненциальное спадание в свободном пространстве от границы метаматериала, поэтому приближенно можно считать, что $h=0$ на образующих клина, и искать собственное решение в виде $h=H(r)\cos(\pi\phi/2\phi_{0})$, где $2\phi_{0}$ угол раскрыва клина и угол $\phi$ отсчитывается от диагонали клина. Введем новую радиальную переменную $\zeta=\log r$ и запишем уравнение (3) с учетом выбранной угловой зависимости

[TABLE]

Первые два члена в (4) дают уравнение осциллятора (напомним, что $\varepsilon_{\phi\phi}<0$), а последнее слагаемое мало при малых $r$ и быстро растет с увеличением расстояния от вершины клина. Критическое значение $r_{max}$, при котором последние два слагаемых становятся величинами одного порядка, можно оценить как

[TABLE]

а при малых $r<r_{max}$ справедливо приближенное выражение

[TABLE]

В соответствии с этой формулой при уменьшении угла раскрыва клина положение максимума поля $H$ сдвигается в сторону больших $r$, что иллюстрируется на рис. 6, где также можно заметить рост числа осцилляций поля. Очевидно, что число осцилляций при уменьшении $r$ стремится к бесконечности, а их характерный масштаб к нулю. Естественные ограничения на минимальный масштаб накладываются физическими причинами, т.е. поглощением в среде, шероховатостью границы метаматериала, его неоднородностью, и др. По-видимому, возбуждение этой собственной моды наблюдалось также в численном эксперименте cui_NanoLett12 , где также указывалось на субволновую локализацию и эффективное поглощение излучения в малом объеме.

V Заключение

В качестве заключения, в работе предложен метод расчета полей, рассеянных на профилированной границе двух однородных анизотропных сред для случая, когда профиль границы является произвольной периодической функцией координаты вдоль границы. Естественное обобщение второго тождества Грина позволило рассмотреть нетривиальный случай, когда одна из сред представляет собой гиперболический метаматериал, характеризующийся условием $\varepsilon_{XX}\varepsilon_{YY}<0$. В этом случае перенормировка координат превращает уравнение гиперболического типа в обычное уравнение Гельмгольца, записанное однако в комплексных пространственных переменных. Метод работает и в поглощающей/активной среде, когда сами компоненты диэлектрического тензора являются комплексными величинами.

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

Найдено приближенное аналитическое выражение для собственной моды, локализованной в гиперболическом метаматериале, который имеет пилообразный профиль границы. Теоретически сколь угодно малый масштаб локализации моды ограничивается диссипацией в среде. Определены параметры (ориентация оси анизотропии, угол раскрыва "пилы"), при которых такое (супер)субволновое решение реализуется.

Работа выполнена при поддержке РФФИ (грант 16-02-00556).

Список литературы

• (1)

A. Poddubny, I. Iorsh, P. Belov, Y. Kivshar, Nat. Photonics 7, 948 (2013).

• (2) D. Lu, Z. Liu, Nat. Commun. 3, 1205 (2012).
• (3) I. V. Iorsh, I. S. Mukhin, I. V. Shadrivov, P. A. Belov, Y. S. Kivshar, Phys. Rev. B 87, 075416 (2013).
• (4) Y. Guo, C. L. Cortes, S. Molesky, and Z. Jacob, Appl. Phys. Lett., 101, 131106 (2012).
• (5)

S. M. Hashemi, I. S. Nefedov, Phys. Rev. B 86, 195411 (2012).

• (6)

I. S. Nefedov, C. A. Valagiannopoulos, S. M. Hashemi, E. I. Nefedov, Sci. Rep. 3, 2662 (2013).

• (7)

I. S. Nefedov, C. A. Valagiannopoulos, L. A. Melnikov, J. Opt. 15, 114003 (2013).

• (8) J. Yao, Z. Liu, Y. Liu, et. al., Science 321, 930 (2008).
• (9) A. J. Hoffman, L. Alekseyev, S. S. Howard, K. J. Franz, et. al., Nat. Mater. 6, 946 (2007).
• (10) Z. Jacob, L. V. Alekseyev, E. Narimanov, Opt. Express 14, 8247 (2006).
• (11)

J. Sun, M. Shalaev, N. Litchinitser, Nat. Commun. 6, 7201 (2015).

• (12)

T. U. Tumkur, L. Gu,J. K. Kitur, E. E. Narimanov, M. A. Noginov, Appl. Phys. Lett. 100, 161103 (2012).

• (13)

T. U. Tumkur, J. K. Kitur, B. Chu, et. al., Appl. Phys. Lett. 101, 091105 (2012)

• (14) Z. Liu, H. Lee, Y. Xiong, C. Sun, X. Zhang, Science 315, 1686 (2007).
• (15)

K. Mantel, D. Bachstein, and U. Peschel, Opt. Lett. 36, 199 (2011).

• (16) Y. Cui, K. H. Fung, J. Xu, et. al., NanoLett. 12, 1443 (2012).
• (17) H. Hu, D. Ji, X. Zeng, K. Liu, Q. Gan, Sci. Rep. 3, 1249 (2013).
• (18)

Z. Jacob, J.-Y. Kim, G. V. Naik, et. al., Appl. Phys. B 100, 215 (2010).

• (19)

P. Shekhar, J. Atkinson, Z. Jacob, Nano Convergence 1, 1 (2014).

• (20) N. A. Zharova, A. A. Zharov, and A. A. Zharov, J. Opt. Soc. Am. B 33, 594 (2016).
• (21) N. A. Zharova, A. A. Zharov , and A. A. Zharov Jr. Advances in Condensed Matter Physics 2018, 4578148, (2018). https://doi.org/10.1155/2018/4578148
• (22) L. Ferrari, C. Wu, D. Lepage, X. Zhang, Z. Liu, Progress in Quantum Electronics 40, 1 (2015).
• (23) www.refractiveindex.info
• (24)

Q. Sun, E. Klaseboer, et.al. R. Soc. open sci. 2, 140520 (2015).

VI Приложение: метод расчета полей при рассеянии плоской волны на периодически модулированной границе гиперболического метаматериала

Метод основан на использовании второго тождества Грина (Green’s second identity) и сводит задачу рассеяния падающего излучения на границе однородной среды к вычислению поверхностных (для двумерного случая, контурных) интегралов от функции Грина и ее производной по нормали к поверхности. Похожая методика применялась в публикации Sun для решения уравнения Гельмгольца в акустике, но в данной работе проведено обобщение, позволяющее рассматривать анизотропную среду гиперболического типа, в которой обычное (эллиптическое) уравнение Гельмгольца фактически становится гиперболическим уравнением.

Второе тождество Грина формулируется следующим образом. Если $\psi$ и $\phi$ обе дважды непрерывно дифференцируемы в области $U\in R^{3}$, тогда

[TABLE]

В рассматриваемом двумерном случае лапласиан $\Delta=\partial^{2}/\partial x^{2}+\partial^{2}/\partial y^{2}$, интегрирование по объему заменяется на интегрирование по $x$, $y$ координатам, а граница области $\partial U$ представляет собой (параметрически заданную) линию $x_{B}(\tau)$, $y_{B}(\tau)$ (см. (1)).

Возьмем в качестве $\phi$ рассеянное магнитное поле $h^{sct}$, а в качестве $\psi$ функцию Грина $G$ двумерной задачи, в которой точка особенности расположена на границе области $U$.

Для излучения TM поляризации уравнение Гельмгольца в свободном пространстве имеет вид

[TABLE]

где $h\equiv h_{z}$ это (единственная) $z$ компонента магнитного поля. Для функции Грина, удовлетворяющей уравнению

[TABLE]

известно аналитическое выражение $G(\mathbf{r},\mathbf{r}_{0})=i/4H^{1}_{0}(k_{0}\rho)$, где $H^{1}_{0}$ функция Ханкеля нулевого порядка и $\rho=\sqrt{(x-x_{0})^{2}+(y-y_{0})^{2}}$.

Нетрудно заметить, что при таком выборе $\phi$ и $\psi$ интеграл в левой части равенства (5) равен нулю, а правая часть задает соотношение между значениями рассеянного магнитного поля $h^{sct}$ и его производной по нормали $\partial h^{sct}/\partial n$ в точках на границе ($x_{B}(\tau_{i})$, $y_{B}(\tau_{i})$), причем аналогичное соотношение можно записать для каждой $i$-той из $N$ граничных точек, используя соответствующие функции Грина, $G_{ij}\equiv G(k_{0}\rho_{ij})$, где $\rho_{ij}=\sqrt{[x_{B}(\tau_{i})-x_{B}(\tau_{j})]^{2}+[y_{B}(\tau_{i})-y_{B}(\tau_{j})]^{2}}$.

Таким образом, для $2N$ неизвестных (значения $h^{sct}$ и $\partial h^{sct}/\partial n$) у нас есть $N$ уравнений, полученных из соотношений (5),

[TABLE]

где $d\tau_{j}=\sqrt{dx_{j}^{2}+dy_{j}^{2}}$ и как обычно предполагается суммирование по повторяющимся индексам. Недостающие $N$ уравнений находятся из условий непрерывности тангенциальных компонент электрического и магнитного полей на границе $x_{B}(\tau_{i})$, $y_{B}(\tau_{i})$. В простейшем случае, когда падающее излучение рассеивается на идеально отражающем металле, эти условия сводятся к равенству нулю нормальной производной полного (рассеянное плюс падающее) магнитного поля, $\partial h^{sct}/\partial n=-\partial h^{inc}/\partial n$.

Если граница $x_{B}(\tau)$, $y_{B}(\tau)$ разделяет свободное пространство и диэлектрик с проницаемостью $\varepsilon_{D}$, то в области диэлектрика волновое число $k=\sqrt{\varepsilon_{D}}k_{0}$, и функция Грина перенормируется как $\tilde{G}=G(k\rho/k_{0})$. Прошедшее магнитное поле $h^{tr}$ и его нормальная производная на границе удовлетворяет уравнению  (5) с этой новой функцией Грина $\tilde{G}$, а условия непрерывности полей записываются как $h^{tr}=h^{sct}+h^{inc}$ and $\partial h^{tr}/\partial n=-\varepsilon_{D}(\partial h^{sct}/\partial n+\partial h^{inc}/\partial n)$ 111 Знак “ - ” здесь учитывает противоположное направление векторов нормали для области свободного пространства и области диэлектрика, см. рис.7. в каждой из $N$ граничных точек. Таким образом, мы имеем $4N$ уравнений для $4N$ неизвестных, что единственным образом определяет решение.

Следует отметить, что при вычислении коэффициентов матричного уравнения (8) (и аналогичного уравнения для поля и его нормальной производной в области диэлектрика) возникает проблема сингулярности функции Грина $G_{ii}$ и ее производной $\partial G_{ii}/\partial n_{i}$. Очевидно, что требуется более точное вычисление диагональных коэффициентов в уравнении (8). Правильный результат дает сдвиг контура интегрирования таким образом, чтобы точка особенности оставалась вне области $U$, как это иллюстрируется на рис.7. Соответственно, интеграл вычисляется вдоль верхней полуокружности для полей в области свободного пространства и вдоль нижней полуокружности для полей в диэлектрике (см. рис. 7).

Задача усложняется в случае, когда рассеяние происходит на границе среды с тензорной проницаемостью $\hat{\varepsilon}$, например, метаматериала. Уравнение Гельмгольца для магнитного поля в такой среде имеет вид

[TABLE]

где $X$ и $Y$ главные оси, в которых диэлектрический тензор диагонален (в общем случае оси $X,Y$ повернуты на некоторый угол $\theta$ относительно осей $x,y$, определяющих ориентацию границы: $X=x\cos\theta+y\sin\theta$, $Y=y\cos\theta-x\sin\theta$), и $\varepsilon_{XX}=\varepsilon_{\|}$, $\varepsilon_{YY}=\varepsilon_{\bot}$.

Нетрудно заметить, что заменой

[TABLE]

уравнение (9) приводится к каноническому виду (6), и функция Грина в этих переменных по-прежнему выражается через функцию Ханкеля нулевого порядка, $\tilde{G}=i/4H^{1}_{0}(k_{0}\tilde{\rho})$, где $\tilde{\rho}=\sqrt{(\xi-\xi_{0})^{2}+(\eta-\eta_{0})^{2}}=\sqrt{\varepsilon_{YY}(X-X_{0})^{2}+\varepsilon_{XX}(Y-Y_{0})^{2}}$. Поэтому естественным для решения задачи представляется следующий алгоритм: в метаматериальной области вводятся по формулам (10) координаты $\xi$, $\eta$, определяются граничные точки $\ \xi_{B}$, $\eta_{B}$, соответствующие реальным граничным точкам $x_{B}$, $y_{B}$ в свободном пространстве, и вектора нормали $\mathbf{\tilde{n}}$ ($\tilde{n}_{\xi}=d\eta/d\tilde{\tau}$, $\tilde{n}_{\eta}=-d\xi/d\tilde{\tau}$, где $d\tilde{\tau}=\sqrt{d\xi^{2}+d\eta^{2}}$) в этом нормированном пространстве. Интегрирование уравнений (5) в $\xi$, $\eta$ пространстве дает, как и раньше, $N$ уравнений, связывающих $h^{tr}$ and $\partial h^{tr}/\partial\tilde{n}=\tilde{n}_{\xi}\partial h^{tr}/\partial\xi+\tilde{n}_{\eta}\partial h^{tr}/\partial\eta$. Заключительным шагом является использование условий непрерывности магнитного и тангенциального электрического полей на границе между областями.

Нас интересует случай, когда диэлектрическая среда является гиперболическим метаматериалом, в котором продольная и поперечная компоненты диэлектрического тензора имеют разные знаки (для определенности будем считать, что $\varepsilon_{XX}>0$, $\varepsilon_{YY}<0$). При этих условиях нормированные координаты $\xi$ оказываются мнимыми, понятие вектора нормали, как вектора единичной длины, теряет смысл. Что более важно, эффективное расстояние до источника в функции Грина $\tilde{\rho}$ стремится к нулю в направлении резонансного конуса (при $(X-X_{0})/(Y-Y_{0})=\sqrt{-\varepsilon_{XX}/\varepsilon_{YY}}$), а поле в этом направлении бесконечно велико. Ситуацию, как всегда, спасает учет затухания в метаматериале, однако при этом комплексными становятся как $\xi$, так и $\eta$, и интегрирование уравнения (5) приходится проводить в комплексном пространстве.

Однако, можно свести вычисление коэффициентов матричного уравнения (8) к интегрированию в обычном пространстве. Действительно, учтем, что

[TABLE]

где $G^{\prime}$ означает производную от функции Грина по аргументу. Аналогично можно записать выражение для величины

[TABLE]

где $E^{tr}_{\tau}$ тангенциальная компонента электрического поля в прошедшей волне в граничных точках.

Таким образом, уравнение (8) модифицируется для случая среды с тензорной диэлектрической проницаемостью следующим образом

[TABLE]

и аналогичное соотношение записывается для рассеянного поля в свободном пространстве

[TABLE]

Условие непрерывности тангенциальных электрического и магнитного полей приводит к равенствам

[TABLE]

Описанная выше постановка задачи непосредственно применима для расчета рассеяния плоской волны на диэлектрическом/метаматериальном объекте конечных размеров. В этом случае граница двух сред является замкнутым контуром $\partial U$ в формуле (5) для метаматериала, а бесконечно удаленная граница свободного пространства не влияет на процесс рассеяния. Если же мы рассматриваем рассеяние падающего излучения на периодическом профиле и хотим ограничиться численными расчетами полей в пределах элементарной ячейки, имея в виду применение в дальнейшем теоремы Блоха, то замкнутый контур $\partial U$ в формуле (5), ограничивающий область свободного пространства, необходимо содержит участки, которые не граничат с областью метаматериала и для которых поэтому неприменимо граничное условие (13). Соответствующие участки имеются также на замкнутом контуре, ограничивающем область метаматериала (см. рис. 7). На этих участках границы следует использовать условия радиационного излучения (условия излучения Зоммерфельда), которые связывают между собой тангенциальные компоненты электрического и магнитного полей (для свободного пространства, магнитное поле и его нормальную производную). Учет этой связи позволяет полностью рассчитать задачу и найти распределение тангенциальных компонент полей на границе $\partial U$ и тем самым коэффициент отражения.

Если же кроме коэффициента отражения нас интересует распределение полей внутри этого замкнутого контура, то для этого необходимо найти амплитуды эффективных источников $V_{i}$, которые мы должны расположить лишь на части контура $\partial U$, разделяющей свободное пространство и метаматериал и которые должны создавать на этой части контура рассчитанное выше распределение полей. Математически задача сводится к решению матричного уравнения

[TABLE]

и аналогично

[TABLE]

где $\mathbf{r}_{i}$, $\mathbf{r}_{j}$ принадлежат граничной части контура $\partial U$. Соответственно поля в части свободного пространства и метаматериала, ограниченных контурами $\partial U$ и $\partial\tilde{U}$ вычисляются как

[TABLE]

[TABLE]

где $h^{inc}(x,y)$ магнитное поле в падающей плоской волне.