Reflection from the periodically modulated interface of a hyperbolic metamaterial
Nina A. Zharova, Alexander A. Zharov, Alexander A. Zharov Jr
aa r X i v : . [ phy s i c s . op ti c s ] J a n Отражение от профилированной границы гиперболического метаматериала
Н. А. Жарова ∗ Институт прикладной физики, РАН, Нижний Новгород, 603950 Россия
А. А. Жаров и А. А. Жаров, мл.
Институт физики микроструктур, РАН, Нижний Новгород, 603950 Россия
Проведено численное исследование особенностей рассеяния падающего излучения на пери-одически модулированной границе гиперболического метаматериала. Найдена оптимальная- пилообразная - форма профиля, обеспечивающая минимум отражения. Показано, что приопределенной ориентации оптических осей метаматериала на пилообразной границе осуществ-ляется субволновая локализация поля, найдено приближенное аналитическое выражение длясобственной моды. Проведено обобщение численного метода расчета полей, основанного наприменении второго тождества Грина.
PACS коды:
ВВЕДЕНИЕ
Большой интерес исследователей к так называемымгиперболическим метаматериалам (ГМ) [1–7] обуслов-лен их необычными и полезными для приложенийэлектромагнитными свойствами. ГМ это анизотроп-ная резонансная среда, которая характеризуется изо-частотной поверхностью гиперболического типа, а неэллипсоидальной, как для обычных сред. В оптикеГМ может быть реализован, например, как металло-диэлектрическая планарная наноструктура, либо какрешетка металлических нанопроволок в диэлектриче-ской матрице.Уникальные свойства ГМ могут найти множествоприложений, среди которых отрицательная рефрак-ция [8, 9], создание субволновых изображений [10–15], локализация света в широком частотном диа-пазоне [16, 17]. Очень высокая плотность фотонныхсостояний в ГМ [18–20] обеспечивает эффективныйизлучательный теплообмен [4–6], что может исполь-зоваться в дизайне наноразмерных устройств нагре-ва/охлаждения.Гиперболический характер метаматериала достига-ется лишь при достаточно высокой объемной доле ме-талла в композите, что также означает большие оми-ческие потери. Поглощение излучения особенно силь-но для мелкомасштабных мод, когда фазовая ско-рость волны направлена вдоль образующей резонанс-ного конуса. Для оптической литографии и созданияизображений с субволновым разрешением омическиепотери нежелательны, но сильное поглощение в ГМможет быть использовано для дизайна поглощающихпокрытий. В работе [21] был предложен дизайн погло-щающего покрытия на основе неоднородного гипер-болического метаматериала. Оказывается, что полноепоглощение в широкой полосе частот можно реали-зовать в плоском тонком слое ГМ при условии, чтоориентация оптической оси плавно меняется от па-раллельной к нормальной (по отношению к интерфей- су). Однако к сожалению в таком тонком слое трудноконтролировать направление оптической оси метама-териала (т.е., например, ориентацию металлическихпроволочных включений).Идеальный поглотитель должен обеспечивать нетолько полное поглощение излучения в возможно бо-лее тонком слое покрытия. Очень важно достигнутьминимального отражения на границе с воздухом. Из-вестный способ снижения коэффициента отраженияэто структурирование поверхности поглощающего ма-териала, то есть профилирование его границы с возду-хом. В данной работе мы изучаем особенности рассе-яния падающего излучения на слое гиперболическогометаматериала с профилированной границей, в част-ности влияние формы, продольного и поперечного ха-рактерного масштаба неоднородности профиля на ко-эффициент отражения.
ПОСТАНОВКА ЗАДАЧИ
В простейшей постановке задача сводится к отра-жению плоской волны TM поляризации от двумерно-го профиля ˆ ε . Метаматериал описывается диэлектри-ческим тензором ˆ ε , и в главных осях ˆ ε = [ ε k ,
0; 0 , ε ⊥ ] .Граница x B ( τ ) , y B ( τ ) , разделяющая свободное про-странство и метаматериал, является периодическойфункцией параметра τ (длина вдоль границы) и ха-рактеризуется масштабом вдоль границы L x , ампли-тудой профиля L y и формой. Масштабирование про-филя по x и y координатам для заданной формы вы-полняется элементарно, но задача выбора формы нетак тривиальна. Оптимальным представляется зада-ние x B ( τ ) , y B ( τ ) математической формулой с пара-метром, изменение которого существенно меняет фор-му границы. В дальнейшем мы использовали функ-цию, задающую x B ( τ ) , y B ( τ ) на длине элементарнойячейки τ = [ − π/ / π ] , в виде x B = π Z cos θdτ − P sin(2 τ ) , y B = π Z sin θdτ, (1)где θ = − . π cos τ . Управляющий параметр P меняетформу профиля от почти пилообразной к почти ко-синусоидальной и до напоминающей меандр. На ри-сунке 1 приведены в качестве иллюстрации кривые y B ( x B ) для набора значений P . P=−1.2 P=−0.2 P=0.8 P=1.8 P=2.8 P=3.8
Рис. 1: (color online) Форма профиля y B ( x B ) периодиче-ской границы (элементарная ячейка) для различных зна-чений параметра P . Перебирая различные значения параметра P , мас-штаба элементарной ячейки вдоль границы L x и ам-плитуды профиля L y , можно найти оптимальное соче-тание, которое обеспечивает слабое отражение падаю-щего излучения от метаматериала. Дополнительнымпараметром оптимизации может служить также уголориентации главных (оптических) осей метаматериа-ла по отношению к осям x , y .Коэффициент отражения от полупространства ГМс профилированной границей может быть найденпри решении уравнения Гельмгольца методом конеч-ных элементов (FEM). Однако попытка использоватьстандартную FEM программу решения уравнений вчастных производных оказалась неудачной. Посколь-ку в этом случае нет прямого доступа к коду, то возни-кают проблемы: появляется отражение излучения отграниц расчетного интервала, приходится рассматри-вать отражение не плоской волны, но широкого гаус-сова пучка из-за того, что во встроенном модуле нетопции задания периодических по x граничных усло-вий. Это все вносит ошибки в расчет коэффициентаотражения. Более удобным оказался поэтому другойметод расчета, основанный на использовании второготождества Грина (Green’s second identity), подробноеописание которого дается в Приложении. РЕЗУЛЬТАТЫ ЧИСЛЕННОГОИССЛЕДОВАНИЯ
Тестирование метода и соответствующей расчетнойпрограммы проводилось для предельных случаев (i) слабой модуляции границы, (ii) одинаковых диэлек-трических характеристик двух сред, (iii) малого ха-рактерного периода модуляции, когда работает при-ближение эффективной среды. Расчетные и теорети-ческие значения коэффициента отражения совпали вэтих случаях с точностью ∼ % (ошибка обусловленадискретизацией вычислений). Кроме того, мы прово-дили сравнение распределения полей, полученного на-шим методом и рассчитанного стандартной програм-мой, которая использует метод конечных элементов.С учетом сказанного выше относительно трудностейприменения стандартной программы, это сравнениепоказало разумное совпадение.Параметры гиперболической среды, использован-ные в численном моделировании были найдены дляметаматериала представляющего собой диэлектриче-скую матрицу с металлическими включениями в фор-ме проволок (wire medium). В приближении эффек-тивной среды компоненты диэлектрического тензора ˆ ε вдоль проволок ε k и поперек них ε ⊥ находятся как[22] ε k = ε M ρ M + ε D (1 − ρ M ) ,ε ⊥ = ε D ε M (1 + ρ M ) + ε D (1 − ρ M ) ε D (1 + ρ M ) + ε M (1 − ρ M ) , (2)где ε M диэлектрическая проницаемость металла, ε D проницаемость диэлектрической матрицы и ρ M объ-емная доля металла в составе композита. Далее мы ис-пользуем для расчетов параметры стекла ( ε D =2.1590)и серебра с объемной долей ρ M = 0 . . Дисперсионнаязависимость проницаемости серебра ε M ( λ ) берется избазы данных [23]. Для вакуумной длины волны излу-чения λ = 0 . µ m величина ε M = − . . i , ирасчеты по формулам (2) дают значения ε k = 4 . . i и ε ⊥ = − . . i .Результаты расчета рассеяния излучения при нор-мальном падении на периодическую границу ГМ с по-мощью метода, описанного в Приложении, иллюстри-руются на рисунке 2, где приведены значения коэф-фициента отражения как функции амплитуды моду-ляции A/λ и нормированного волнового числа моду-ляции κ /k = L x /λ для двух различных типов про-филя, пилообразного (с параметром формы P = − . )и косинусоидального (с параметром формы P = 0 . ).Если для пилообразного профиля рост амплитуды мо-дуляции приводит к снижению отражения по сравне-нию со случаем плоской границы, то максимальныйкоэффициент отражения от косинусоидального про-филя превышает его почти в два раза. Такое нетриви-альное поведение коэффициента отражения обуслов-лено тензорным и, главным образом, гиперболиче-ским характером среды.На рисунке 3 приведены рассчитанные распределе-ния поля для трех различных типов профиля грани-цы гиперболической среды. Коэффициент отражения (cid:79) (cid:78) / k (cid:79) a b Рис. 2: (color online) Коэффициент отражения от профи-лированной границы с параметром формы P = − . (a)и P = 0 . (b) как функция амплитуды модуляции A/λ инормированного волнового числа модуляции κ /k = L x /λ . плоской волны h inc = exp( − ik y ) оказывается рав-ным R = ε = ε k , но максимальное отражение наблюдается отмеандроподобного профиля (панель(с)), и коэффици-ент отражения в этом случае заметно больше, чем отплоской границы. Рис. 3: (color online) Распределение модуля магнитного по-ля | h ( x, y ) | (показан один период структры) при нормаль-ном падении плоской волны h inc = exp ( − ik y ) на периоди-чески модулированную границу гиперболического метама-териала. Амплитуда модуляции A = 0 . λ , период структу-ры в x направлении λ/ . , параметр P , контролирующийформу модуляции (см. формулу (1)) P = − . , 0.5 и 3.5 со-ответственно для панелей (a), (b) и (c). Рассчитанный дляэтих трех вариантов коэффициент отражения R = Коэффициент отражения существенно зависит оториентации осей анизотропии метаматериала по от-ношению к границе. При повороте осей на 90 o , ко- гда компонента тензора проницаемости вдоль грани-цы отрицательна, ε xx = ε ⊥ = − . . i , ко-эффициент отражения от плоской границы близок кединице, и профилирование границы уменьшает от-ражение. На рисунке 4 приведена зависимость коэф-фициента отражения от профилированной границы спараметром формы P = − . (a) и P = 0 . (b) какфункция амплитуды модуляции A/λ и нормирован-ного волнового числа модуляции κ /k = L x /λ . Су-щественное уменьшение отражения наблюдается приусловии роста как амплитуды модуляции A , так и вол-нового числа модуляции κ для P = − . . Эта же тен-денция, хотя и с гораздо меньшим эффектом, просле-живается для профиля другой формы с параметром P = 0 . (форма профиля изображена на рис. 3 (b)). (cid:79) (cid:78) / k (cid:79) a b Рис. 4: (color online) Коэффициент отражения от профи-лированной границы с параметром формы P = − . (a)и P = 0 . (b) как функция амплитуды модуляции A/λ инормированного волнового числа модуляции κ /k = L x /λ .Ориентация осей анизотропии повернута на 90 o по срав-нению со случаем изображенном на рисунке 2 Следует отметить, что основным параметром, опре-деляющим отражение (по крайней мере, при доста-точно сильной модуляции границы) является отно-шение поперечного (амплитуды модуляции A ) и про-дольного ( L x = 2 π/ κ ) масштабов модуляции, то естьпараметр g = A κ / π . На рисунке 5 приведена за-висимость R ( g ) , где красными кружками (чернымизвездочками) отмечены данные для пилообразного, P = − . (косинусоидального, P = 0 . ) профиля, ималая степень дисперсии коэффициента отраженияпри g > . указывает на то, что одинаковый эффектдостигается как при сильной (относительно длинно-волновой), так и слабой (но мелкомасштабной) моду-ляции. СУБВОЛНОВАЯ СОБСТВЕННАЯ МОДА
Интересный эффект наблюдается для такой ориен-тации оптических осей ГМ при рассеянии падающегоизлучения на пилообразом профиле. При достаточ- κ /2 π Рис. 5: (color online) Коэффициент отражения от профили-рованной границы с параметром формы P = − . (крас-ные кружки) и P = 0 . (черные звездочки) как функцияпараметра A κ / π . Ориентация осей анизотропии повер-нута на 90 o по сравнению со случаем изображенном нарисунке 2 но узких зубцах “пилы” поле локализуется вблизи еевершины. Характерные распределения поля приведе-ны на рисунке 6. Эти локализованные распределения Рис. 6: (color online) Распределение модуля магнитного по-ля h ( x, y ) (показан один период структры) при нормаль-ном падении плоской волны h inc = exp ( − ik y ) на пери-одически модулированную границу гиперболического ме-таматериала. Амплитуда модуляции A = 0 . λ , параметр P , контролирующий форму модуляции (см. формулу (1)) P = − . , период структуры в x направлении λ/ κ , где κ /k =
5, 7, и 9 соответственно для панелей (a), (b) и (c).Ориентация осей анизотропии повернута на 90 o по сравне-нию со случаем изображенном на рисунке 3. Коэффициентотражения для этих вариантов 1%, 4.3% и 0.9% соответ-ственно. имеют субволновой характерный масштаб, и анализпоказывает, что теоретически масштаб локализацииможет быть сколь угодно малым.Действительно, в этом случае задача имеет прибли-женное аналитическое решение. Для анализа перей-дем в полярную систему координат ( r , φ ), с центромв вершине клина (зуба пилы) и будем считать, что ди-электрический тензор также имеет полярную симмет-рию, так что диагональные компоненты ε φφ = ε ⊥ < и ε rr = ε k > . Уравнение Гельмгольца для магнитно-го поля в этом случае имеет вид ε − φφ r ∂∂r r ∂h∂r + ε − rr r ∂ h∂φ + k h = 0 . (3) Очевидно, что субволновая локализация поля пред-полагает его экспоненциальное спадание в свобод-ном пространстве от границы метаматериала, поэто-му приближенно можно считать, что h = 0 на обра-зующих клина, и искать собственное решение в виде h = H ( r ) cos( πφ/ φ ) , где φ угол раскрыва клинаи угол φ отсчитывается от диагонали клина. Введемновую радиальную переменную ζ = log r и запишемуравнение (3) с учетом выбранной угловой зависимо-сти d Hdζ − ε φφ ε rr (cid:18) π φ (cid:19) H + k ε φφ e ζ H = 0 . (4)Первые два члена в (4) дают уравнение осциллятора(напомним, что ε φφ < ), а последнее слагаемое малопри малых r и быстро растет с увеличением рассто-яния от вершины клина. Критическое значение r max ,при котором последние два слагаемых становятся ве-личинами одного порядка, можно оценить как k r max ≈ ε − / rr π φ , а при малых r < r max справедливо приближенное вы-ражение H ≈ sin (cid:18)q − ε φφ /ε rr π φ log r (cid:19) . В соответствии с этой формулой при уменьшении уг-ла раскрыва клина положение максимума поля H сдвигается в сторону больших r , что иллюстрируетсяна рис. 6, где также можно заметить рост числа ос-цилляций поля. Очевидно, что число осцилляций приуменьшении r стремится к бесконечности, а их харак-терный масштаб к нулю. Естественные ограниченияна минимальный масштаб накладываются физически-ми причинами, т.е. поглощением в среде, шероховато-стью границы метаматериала, его неоднородностью, идр. По-видимому, возбуждение этой собственной мо-ды наблюдалось также в численном эксперименте [16],где также указывалось на субволновую локализациюи эффективное поглощение излучения в малом объе-ме. ЗАКЛЮЧЕНИЕ
В качестве заключения, в работе предложен методрасчета полей, рассеянных на профилированной гра-нице двух однородных анизотропных сред для слу-чая, когда профиль границы является произвольнойпериодической функцией координаты вдоль грани-цы. Естественное обобщение второго тождества Гринапозволило рассмотреть нетривиальный случай, когдаодна из сред представляет собой гиперболический ме-таматериал, характеризующийся условием ε XX ε Y Y < . В этом случае перенормировка координат пре-вращает уравнение гиперболического типа в обыч-ное уравнение Гельмгольца, записанное однако в ком-плексных пространственных переменных. Метод ра-ботает и в поглощающей/активной среде, когда самикомпоненты диэлектрического тензора являются ком-плексными величинами.C помощью разработанного метода численно иссле-дованы особенности рассеяния падающего излученияна слое гиперболического метаматериала с периоди-чески профилированной границей, в частности, влия-ние формы, продольного и поперечного характерногомасштаба неоднородности профиля на коэффициентотражения. Численное моделирование показало, чтодостаточно мелкомасштабное профилирование грани-цы приводит к уменьшению отражения по сравнениюс плоской границей практически для всех рассмот-ренных форм профиля (при изменении управляющегопараметра форма профиля меняется от почти пилооб-разной к почти косинусоидальной и до напоминающеймеандр). Однако чем больше форма профиля откло-няется от пилообразной, тем чувствительнее оказыва-ется коэффициент отражения к изменению амплиту-ды и характерного продольного масштаба модуляции.Найдено приближенное аналитическое выражениедля собственной моды, локализованной в гиперболи-ческом метаматериале, который имеет пилообразныйпрофиль границы. Теоретически сколь угодно малыймасштаб локализации моды ограничивается диссипа-цией в среде. Определены параметры (ориентация осианизотропии, угол раскрыва "пилы"), при которыхтакое (супер)субволновое решение реализуется.Работа выполнена при поддержке РФФИ (грант 16-02-00556). ∗ Electronic address: [email protected][1] A. Poddubny, I. Iorsh, P. Belov, Y. Kivshar, Nat.Photonics , 948 (2013).[2] D. Lu, Z. Liu, Nat. Commun. , 1205 (2012).[3] I. V. Iorsh, I. S. Mukhin, I. V. Shadrivov, P. A. Belov,Y. S. Kivshar, Phys. Rev. B , 075416 (2013).[4] Y. Guo, C. L. Cortes, S. Molesky, and Z. Jacob, Appl.Phys. Lett., , 131106 (2012).[5] S. M. Hashemi, I. S. Nefedov, Phys. Rev. B , 195411(2012).[6] I. S. Nefedov, C. A. Valagiannopoulos, S. M. Hashemi,E. I. Nefedov, Sci. Rep. , 2662 (2013).[7] I. S. Nefedov, C. A. Valagiannopoulos, L. A. Melnikov,J. Opt. , 114003 (2013).[8] J. Yao, Z. Liu, Y. Liu, et. al., Science , 930 (2008).[9] A. J. Hoffman, L. Alekseyev, S. S. Howard, K. J. Franz,et. al., Nat. Mater. , 946 (2007).[10] Z. Jacob, L. V. Alekseyev, E. Narimanov, Opt. Express , 8247 (2006).[11] J. Sun, M. Shalaev, N. Litchinitser, Nat. Commun. ,7201 (2015). [12] T. U. Tumkur, L. Gu,J. K. Kitur, E. E. Narimanov,M. A. Noginov, Appl. Phys. Lett. , 161103 (2012).[13] T. U. Tumkur, J. K. Kitur, B. Chu, et. al., Appl. Phys.Lett. , 091105 (2012)[14] Z. Liu, H. Lee, Y. Xiong, C. Sun, X. Zhang, Science ,1686 (2007).[15] K. Mantel, D. Bachstein, and U. Peschel, Opt. Lett. ,199 (2011).[16] Y. Cui, K. H. Fung, J. Xu, et. al., NanoLett. , 1443(2012).[17] H. Hu, D. Ji, X. Zeng, K. Liu, Q. Gan, Sci. Rep. , 1249(2013).[18] Z. Jacob, J.-Y. Kim, G. V. Naik, et. al., Appl. Phys. B , 215 (2010).[19] P. Shekhar, J. Atkinson, Z. Jacob, Nano Convergence ,1 (2014).[20] N. A. Zharova, A. A. Zharov, and A. A. Zharov, J. Opt.Soc. Am. B , 594 (2016).[21] N. A. Zharova, A. A. Zharov , and A. A. Zharov Jr.Advances in Condensed Matter Physics , 4578148,(2018). https://doi.org/10.1155/2018/4578148[22] L. Ferrari, C. Wu, D. Lepage, X. Zhang, Z. Liu, Progressin Quantum Electronics , 140520(2015). ПРИЛОЖЕНИЕ: МЕТОД РАСЧЕТА ПОЛЕЙПРИ РАССЕЯНИИ ПЛОСКОЙ ВОЛНЫ НАПЕРИОДИЧЕСКИ МОДУЛИРОВАННОЙГРАНИЦЕ ГИПЕРБОЛИЧЕСКОГОМЕТАМАТЕРИАЛА
Метод основан на использовании второго тождестваГрина (Green’s second identity) и сводит задачу рас-сеяния падающего излучения на границе однороднойсреды к вычислению поверхностных (для двумерно-го случая, контурных) интегралов от функции Гринаи ее производной по нормали к поверхности. Похожаяметодика применялась в публикации [24] для решенияуравнения Гельмгольца в акустике, но в данной рабо-те проведено обобщение, позволяющее рассматриватьанизотропную среду гиперболического типа, в кото-рой обычное (эллиптическое) уравнение Гельмгольцафактически становится гиперболическим уравнением.Второе тождество Грина формулируется следую-щим образом. Если ψ и φ обе дважды непрерывнодифференцируемы в области U ∈ R , тогда Z U ( ψ ∆ φ − φ ∆ ψ ) dV = I ∂U (cid:18) ψ ∂φ∂ n − φ ∂ψ∂ n (cid:19) dS. (5)В рассматриваемом двумерном случае лапласиан ∆ = ∂ /∂x + ∂ /∂y , интегрирование по объему заменяет-ся на интегрирование по x , y координатам, а границаобласти ∂U представляет собой (параметрически за-данную) линию x B ( τ ) , y B ( τ ) (см. (1)).Возьмем в качестве φ рассеянное магнитное поле h sct , а в качестве ψ функцию Грина G двумерной за-дачи, в которой точка особенности расположена награнице области U .Для излучения TM поляризации уравнение Гельм-гольца в свободном пространстве имеет вид ∆ h + k h = 0 , (6)где h ≡ h z это (единственная) z компонента магнитно-го поля. Для функции Грина, удовлетворяющей урав-нению ∆ G + k G = δ ( x − x ) δ ( y − y ) ≡ δ ( r − r ) , (7)известно аналитическое выражение G ( r , r ) = i/ H ( k ρ ) , где H функция Ханкеля нулевогопорядка и ρ = p ( x − x ) + ( y − y ) .Нетрудно заметить, что при таком выборе φ и ψ ин-теграл в левой части равенства (5) равен нулю, а пра-вая часть задает соотношение между значениями рас-сеянного магнитного поля h sct и его производной понормали ∂h sct /∂n в точках на границе ( x B ( τ i ) , y B ( τ i ) ),причем аналогичное соотношение можно записать длякаждой i -той из N граничных точек, используя со-ответствующие функции Грина, G ij ≡ G ( k ρ ij ) , где ρ ij = p [ x B ( τ i ) − x B ( τ j )] + [ y B ( τ i ) − y B ( τ j )] .Таким образом, для N неизвестных (значения h sct и ∂h sct /∂n ) у нас есть N уравнений, полученных изсоотношений (5), ∂G ij ∂n j h sctj dτ j − G ij ∂h sctj ∂n j dτ j = 0 , (8)где dτ j = q dx j + dy j и как обычно предполагаетсясуммирование по повторяющимся индексам. Недоста-ющие N уравнений находятся из условий непрерывно-сти тангенциальных компонент электрического и маг-нитного полей на границе x B ( τ i ) , y B ( τ i ) . В простей-шем случае, когда падающее излучение рассеиваетсяна идеально отражающем металле, эти условия сво-дятся к равенству нулю нормальной производной пол-ного (рассеянное плюс падающее) магнитного поля, ∂h sct /∂n = − ∂h inc /∂n .Если граница x B ( τ ) , y B ( τ ) разделяет свободное про-странство и диэлектрик с проницаемостью ε D , то вобласти диэлектрика волновое число k = √ ε D k , ифункция Грина перенормируется как ˜ G = G ( kρ/k ) .Прошедшее магнитное поле h tr и его нормальная про-изводная на границе удовлетворяет уравнению (5) сэтой новой функцией Грина ˜ G , а условия непрерыв-ности полей записываются как h tr = h sct + h inc and ∂h tr /∂n = − ε D ( ∂h sct /∂n + ∂h inc /∂n ) [ ? ] в каждойиз N граничных точек. Таким образом, мы имеем N уравнений для N неизвестных, что единственным об-разом определяет решение.Следует отметить, что при вычислении коэффи-циентов матричного уравнения (8) (и аналогичного −0.2 −0.1 0 0.1 0.2−0.2−0.15−0.1−0.0500.050.10.150.2 aira −0.4 −0.3 −0.2 −0.1 0−0.2−0.15−0.1−0.0500.050.10.150.2 MMb
Рис. 7: (color online) Геометрия задачи: (a) замкнутыйконтур, используемый для расчета рассеянных полей всвободном пространстве, показано направление нормалей,стрелки указывают направление обхода контура. На синейчасти контура применяются условия радиационного излу-чения (см. текст); (b) то же для области метаматериала. уравнения для поля и его нормальной производной вобласти диэлектрика) возникает проблема сингуляр-ности функции Грина G ii и ее производной ∂G ii /∂n i .Очевидно, что требуется более точное вычислениедиагональных коэффициентов в уравнении (8). Пра-вильный результат дает сдвиг контура интегрирова-ния таким образом, чтобы точка особенности остава-лась вне области U , как это иллюстрируется на рис.7.Соответственно, интеграл вычисляется вдоль верхнейполуокружности для полей в области свободного про-странства и вдоль нижней полуокружности для полейв диэлектрике (см. рис. 7).Задача усложняется в случае, когда рассеяние про-исходит на границе среды с тензорной проницаемо-стью ˆ ε , например, метаматериала. Уравнение Гельм-гольца для магнитного поля в такой среде имеет вид ε Y Y ∂ h tr ∂X + 1 ε XX ∂ h tr ∂Y + k h tr = 0 , (9)где X и Y главные оси, в которых диэлектрическийтензор диагонален (в общем случае оси X, Y поверну-ты на некоторый угол θ относительно осей x, y , опре-деляющих ориентацию границы: X = x cos θ + y sin θ , Y = y cos θ − x sin θ ), и ε XX = ε k , ε Y Y = ε ⊥ .Нетрудно заметить, что заменой ξ = √ ε Y Y
X, η = √ ε XX Y (10)уравнение (9) приводится к каноническому виду (6),и функция Грина в этих переменных по-прежнемувыражается через функцию Ханкеля нулевого поряд-ка, ˜ G = i/ H ( k ˜ ρ ) , где ˜ ρ = p ( ξ − ξ ) + ( η − η ) = p ε Y Y ( X − X ) + ε XX ( Y − Y ) . Поэтому естествен-ным для решения задачи представляется следующийалгоритм: в метаматериальной области вводятся поформулам (10) координаты ξ , η , определяются гра-ничные точки ξ B , η B , соответствующие реальнымграничным точкам x B , y B в свободном пространстве,и вектора нормали ˜n ( ˜ n ξ = dη/d ˜ τ , ˜ n η = − dξ/d ˜ τ , где d ˜ τ = p dξ + dη ) в этом нормированном простран-стве. Интегрирование уравнений (5) в ξ , η простран-стве дает, как и раньше, N уравнений, связывающих h tr and ∂h tr /∂ ˜ n = ˜ n ξ ∂h tr /∂ξ + ˜ n η ∂h tr /∂η . Заклю-чительным шагом является использование условийнепрерывности магнитного и тангенциального элек-трического полей на границе между областями.Нас интересует случай, когда диэлектрическая сре-да является гиперболическим метаматериалом, в ко-тором продольная и поперечная компоненты диэлек-трического тензора имеют разные знаки (для опреде-ленности будем считать, что ε XX > , ε Y Y < ). Приэтих условиях нормированные координаты ξ оказыва-ются мнимыми, понятие вектора нормали, как векто-ра единичной длины, теряет смысл. Что более важ-но, эффективное расстояние до источника в функцииГрина ˜ ρ стремится к нулю в направлении резонансно-го конуса (при ( X − X ) / ( Y − Y ) = p − ε XX /ε Y Y ), аполе в этом направлении бесконечно велико. Ситуа-цию, как всегда, спасает учет затухания в метамате-риале, однако при этом комплексными становятся как ξ , так и η , и интегрирование уравнения (5) приходитсяпроводить в комплексном пространстве.Однако, можно свести вычисление коэффициентовматричного уравнения (8) к интегрированию в обыч-ном пространстве. Действительно, учтем, что ∂ ˜ G ij ∂ ˜ n j d ˜ τ j = ˜ G ′ ij k ˜ ρ ij [( ξ j − ξ i ) dη j − ( η j − η i ) dξ j ] = √ ε XX ε Y Y ˜ G ′ ij k ˜ ρ ij [( x j − x i ) dy j − ( y j − y i ) dx j ] = √ ε k ε ⊥ ˜ G ′ ij k ˜ ρ ij ( r j − r i ) n j dτ j , где G ′ означает производную от функции Грина по ар-гументу. Аналогично можно записать выражение для величины ∂h trj ∂ ˜ n j d ˜ τ j = ∂h trj ∂ξ dη j − ∂h trj ∂η dξ j = √ ε k ε ⊥ (cid:18) ∂h trj ∂X dY j ε Y Y − ∂h trj ∂Y dX j ε XX (cid:19) = √ ε XX ε Y Y ik (cid:0) E trX dX j + E trY dY j (cid:1) = √ ε k ε ⊥ ik E trτ dτ j , где E trτ тангенциальная компонента электрическогополя в прошедшей волне в граничных точках.Таким образом, уравнение (8) модифицируется дляслучая среды с тензорной диэлектрической проница-емостью следующим образом X j " ˜ G ′ ij k ˜ ρ ij ( r j − r i ) n j h trj dτ j − ik ˜ G ij ( E trτ ) j dτ j = 0 , (11)и аналогичное соотношение записывается для рассе-янного поля в свободном пространстве X j (cid:20) G ′ ij k ρ ij ( r j − r i ) n j h sctj dτ j − ik G ij ( E sctτ ) j dτ j (cid:21) = 0 . (12)Условие непрерывности тангенциальных электриче-ского и магнитного полей приводит к равенствам h trj = h sctj + h incj ; ( E trτ ) j = ( E sctτ ) j + ( E incτ ) j . (13)Описанная выше постановка задачи непосредствен-но применима для расчета рассеяния плоской волнына диэлектрическом/метаматериальном объекте ко-нечных размеров. В этом случае граница двух средявляется замкнутым контуром ∂U в формуле (5) дляметаматериала, а бесконечно удаленная граница сво-бодного пространства не влияет на процесс рассея-ния. Если же мы рассматриваем рассеяние падаю-щего излучения на периодическом профиле и хотимограничиться численными расчетами полей в преде-лах элементарной ячейки, имея в виду применениев дальнейшем теоремы Блоха, то замкнутый контур ∂U в формуле (5), ограничивающий область свободно-го пространства, необходимо содержит участки, кото-рые не граничат с областью метаматериала и для ко-торых поэтому неприменимо граничное условие (13).Соответствующие участки имеются также на замкну-том контуре, ограничивающем область метаматериа-ла (см. рис. 7). На этих участках границы следует ис-пользовать условия радиационного излучения (усло-вия излучения Зоммерфельда), которые связываютмежду собой тангенциальные компоненты электриче-ского и магнитного полей (для свободного простран-ства, магнитное поле и его нормальную производную).Учет этой связи позволяет полностью рассчитать за-дачу и найти распределение тангенциальных компо-нент полей на границе ∂U и тем самым коэффициентотражения.Если же кроме коэффициента отражения нас инте-ресует распределение полей внутри этого замкнуто-го контура, то для этого необходимо найти амплиту-ды эффективных источников V i , которые мы должнырасположить лишь на части контура ∂U , разделяю-щей свободное пространство и метаматериал и кото-рые должны создавать на этой части контура рассчи-танное выше распределение полей. Математически за-дача сводится к решению матричного уравнения V i G ij = h sctj , и аналогично ˜ V i ˜ G ij = h trj , где r i , r j принадлежат граничной части контура ∂U .Соответственно поля в части свободного пространстваи метаматериала, ограниченных контурами ∂U и ∂ ˜ U вычисляются как h ( x, y ) = V i G ( k p ( x − x i ) + ( y − y i ) ) + h inc ( x, y ) , ˜ h ( X, Y ) = ˜ V i G ( k p ε Y Y ( X − X i ) + ε XX ( Y − Y i ) ) , где h inc ( x, y ))