ЖЭТФ, 2021, том 159, вып. 2, стр. 195-215

ВОЗМОЖНЫ ЛИ ИЗОТРОПНЫЕ МЕТАМАТЕРИАЛЫ

И МЕТАМАТЕРИАЛЫ С ОТРИЦАТЕЛЬНЫМИ

ПРОНИЦАЕМОСТЯМИ ε И μ?

М. В. Давидович^*

Национальный исследовательский Саратовский государственный университет им. Н. Г. Чернышевского

410012, Саратов, Россия

Поступила в редакцию 12 февраля 2020 г.,

после переработки 22 июля 2020 г.

Принята к публикации 26 июля 2020 г.

Предпринята попытка ответа на поставленный в заглавии вопрос. Рассмотрены периодические метамате-

риалы — кубические фотонные кристаллы с диэлектрическими включениями, в том числе с металличес-

кими свойствами, методы введения материальных уравнений (гомогенизации), возможность описания

метаматериалов двумя скалярными проницаемостями (или двумя скалярными параметрами), в том чис-

ле с отрицательными значениями их действительных частей. Для фотонных кристаллов подтвержден вы-

вод Ландау об отсутствии высокочастотных (оптических) магнитных свойств. В частности, невозможны

отрицательные значения магнитной проницаемости или каких-либо ее компонент. Приведены конфигура-

ции метаматериалов, обладающие почти изотропными свойствами при пренебрежении пространственной

дисперсией.

DOI: 10.31857/S0044451021020012

а аналогией для смесей служат аморфные диэлек-

трики. Часто метаматериалы понимаются в узком

1. ВВЕДЕНИЕ

смысле как периодические ИС, содержащие метал-

Под метаматериалами в широком смысле будем

лические МА или частицы, и даже в более узком

понимать любые искусственные среды (ИС), кото-

смысле как ИС с «одновременно отрицательными

рые благодаря развитию современных технологий,

диэлектрической проницаемостью (ДП) ε и магнит-

включая нанотехнологии, в последние два десятиле-

ной проницаемостью (МП) μ» [2, 4]. Также их ас-

тия усиленно создаются и широко теоретически ис-

социируют с ИС, в которых распространяются объ-

следуются. Их исследование началось более ста лет

емные обратные волны (ОВ) [4]. В последнее время

назад, и уже к середине прошлого столетия были

возрос интерес к метаповерхностям, вдоль которых

получены важные результаты по искусственным ди-

возможны ОВ, при этом требование одновременной

электрикам и структурам (см., например, обзор [1]).

отрицательности ДП и МП (или каких-либо их ком-

Создаваемые ИС могут быть периодическими и хао-

понент) не является обязательным [5]. Сразу огово-

тическими, т. е. с периодическим и хаотическим рас-

римся, что ни ε, ни μ отрицательными быть не мо-

положением частиц или метаатомов (МА) в осно-

гут, поскольку это комплексные величины, тем бо-

ве — обычно некой однородной и изотропной ди-

лее что потери в металлических МА существенные,

электрической среде. Другой тип ИС создают и ис-

а рассматриваемые структуры с ОВ резонансные.

следуют на основе теории смесей, формул смеше-

Кроме того, ОВ существуют в любом фотонном кри-

ния, теории перколяции, метода компактных групп

сталле, который строго описывается только тензор-

и ряда других подходов [2, 3]. Здесь мы такие мета-

ной ДП с учетом пространственной дисперсии [6-8].

материалы не рассматриваем, а будем исследовать

Поэтому много работ было посвящено получению и

только периодические ИС, которые также называют

исследованию полностью диэлектрических ИС с од-

фотонными кристаллами. Соответственно имеет ме-

новременно электрическими и магнитными откли-

сто аналогия между ними и кристаллами в оптике,

ками. Ответ на вопрос о том, возможно ли получить

отрицательные эффективные значения ε и μ (а тем

^* E-mail: davidovichmv@info.sgu.ru

195

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

более их одновременно отрицательные значения) в

ционными потерями можно было пренебречь. Имен-

полностью диэлектрических ИС, является одной из

но в этих приближениях мы и будем рассматривать

целей данной работы и будет дан ниже. Обычно

задачу. Например, тонкая пленка в 100 нм с пятью

«одновременность» понимается как значения на од-

периодами и двумя разными слоями по 10 нм в пе-

ной фиксированной частоте. Эта идея заманчива,

риоде не является 1D-кристаллом. Она более соот-

поскольку тангенс угла потерь хороших диэлектри-

ветствует оптическому фильтру [11,12]. Чтобы воз-

ков в широком диапазоне от микроволн до оптики

никло подобие запрещенной зоны, пленка должна

может быть порядка и меньше 10-4, при этом не

иметь толщину более 800 нм, когда число перио-

требуется криогенных температур. Однако эти па-

дов больше или равно 40 [12]. Фотонный кристалл

раметры зависят и от волнового вектора k, т. е. на-

(как и обычный) описывается тензором ДП ε, а ни-

до понимать одновременность как принадлежность

как не отрицательным скаляром. Однако имеется

к фиксированной точке поверхности изочастот (дис-

весьма много работ, в которых периодические ИС

персионной поверхности в k-пространстве). Соот-

описываются с помощью ДП ε < 0 и МП μ < 0.

ветственно переход от одной области дисперсион-

Часто такие параметры просто вводят и анализи-

ной поверхности или ветви к другой может сопро-

руются какие-то свойства, например, поведение лу-

вождаться переходом от прямых волн к обратным.

чей. Тогда волновая векторная электродинамика за-

Снизить потери на два порядка можно, используя

меняется лучевой оптикой. В лучшем случае такие

гелиевые температуры. Такие периодические ИС с

проницаемости вводят в уравнения Максвелла. Это

металлическими МА уже можно рассматривать как

связано с тем, что строгие решения весьма сложны

фотонные кристаллы наряду с диэлектрическими

и в известных нам работах не применялись (за ис-

ИС. В подобных электромагнитных кристаллах по-

ключением использования коммерческих программ-

глощением фотонов можно пренебречь, и тогда по-

ных пакетов). Ниже приведены подходы к строго-

ведение последних аналогично поведению электро-

му численно-аналитическому решению подобных за-

нов в обычных кристаллах. Так же как электроны и

дач.

дырки в твердотельном кристалле, это квазичасти-

Известно, что внутреннее поле играет огромную

цы — квазифотоны или поляритоны, обладающие

роль при описании ДП и МП в природных веще-

дисперсией и позволяющие эффективно управлять

ствах. Обычно ДП и МП вводят путем усреднения

светом, чем и обусловлен большой интерес к таким

полей, индукций и поляризаций микроскопических

структурам нанофотоники. Квазифотоны обладают

уравнений Максвелла - Лоренца по физически бес-

дисперсией, причем для них возможны почти запре-

конечно малому объему [6-8, 13]. Для кристаллов

щенные зоны. Слово «почти» можно опустить, если

и смесей такой процесс называется гомогенизацией.

пренебречь диссипацией и рассматривать бесконеч-

Для смесей усреднение проводится по малому объ-

ный периодический фотонный кристалл.

ему с линейным размером, существенно меньшим

В электродинамике сплошных сред известно

длины волны, но в котором много МА. Для фотон-

утверждение, что с ростом частоты МП μ все бо-

ного кристалла будем делать усреднение по ячейке

лее теряет свой смысл, а в области оптических час-

периодичности. Требовать ее малость по сравнению

тот μ ≈ 1, и рассматривать магнитный отклик бес-

с длиной волны не обязательно, но для локальнос-

смысленно [9] §79, §103. Это утверждение за послед-

ти материальных уравнений необходимо. Низкочас-

ние двадцать лет многократно подвергалось сомне-

тотная гомогенизация означает k = 0, т. е. прене-

нию и опровержению (см., например, [10]). Однако

брежение пространственной дисперсией (ПД). Для

вывод о нецелесообразности высокочастотного опи-

рассчитанной зонной структуры гомогенизация за-

сания природных материалов и метаматериалов с

висит от точки на дисперсионной поверхности. Упо-

помощью МП верен. В данной работе этот вывод,

мянутое усреднение изначально зависит и от модели

полученный в [9] в общем случае, подтверждается

среды [1, 4-8, 14]. Так, можно вовсе не вводить МП,

конкретно для фотонных кристаллов с привлечени-

т. е. не использовать симметричный подход, а опре-

ем классической электродинамики и гомогенизации.

делять только тензор ε. Для определения его шести

Его можно распространить на взаимодействие из-

компонент (в общем случае комплексных) вполне

лучения с периодически расположенными МА как

достаточно шести скалярных уравнений Максвелла,

квантовыми точками. Также рассмотрен вопрос о

что и определяет его преимущества. Симметричный

создании изотропных метаматериалов. Важным в

же подход на основе ε и μ требует введения допол-

этом случае является размер кристалла или число

нительных условий, обычно связанных с симметри-

периодов. Их должно быть так много, чтобы радиа-

ей [15,16].

196

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

2. ИСТОРИЯ ВОПРОСА

В первых работах по гомогенизации [17-19] при

задании модели среды в виде скаляров ε и μ бы-

ли получены магнитные свойства с использовани-

ем только диэлектрических МА. Так, Левин [17, 18]

рассматривал отражение нормально падающей вол-

ны от кубической решетки сферических магнитоди-

электрических МА. Поляризация волны была фик-

сирована вдоль одной из осей. Использовался ряд

приближений. Для металлических МА был получен

диамагнетизм. Для диэлектрических МА с большой

ДП получена широкая вариация значений μ с утвер-

ждением, что эта величина может лежать почти в

любой точке комплексной плоскости. Это не совсем

так. На рис. 1 приведены результаты вычислений по

формулам из [17,18], из которых видно, что μ^′e > 0,

μ^′′e > 0 в области ε′2 < 5000. При более высоких зна-

чениях ДП и достаточно больших шариках возмож-

ны резонансы с малой зоной отрицательной ДП. Од-

Рис. 1. Re (εe) (кривые 1-4) и Re (μe) (5-8) в зависимости

нако веществ с такими большими линейными изо-

от Re (ε2) согласно [12,13] при μ1 = μ2 = 1, k0d√ε1 = 0.1:

тропными ДП нет, а сами резонансы обусловлены

ε1

= 3 (1-3, 5-7), 1 (4, 8). Кривые 1, 5 соответствуют

(

a/d = 0.1; кривые 2, 3 — a/d = 0.2; кривые 4, 7, 8 —

введенной в [18] функцией F

k₀a√ε2μ2) и являют-

a/d = 0.3

ся нефизическими. При больших ДП шариков ДП ε^′e

насыщается, что также говорит об ограниченности

модели. ДП основы в расчетах равнялась 3-i0.0003,

√

при этом получено μ^′′e/μ^′e ∼ 10-4. ДП шариков из-

μ_e/ε_e и n =

√ε_eμ_e. Два условия определяют две

менялась от 1 до 10000, при этом тангенс угла по-

эффективные величины ε_e и μ_e. Эти величины раз-

терь брался равным 10-4. Сам автор отметил, что

личны по разным направлениям. При этом даже

результаты пригодны только для макроскопическо-

если слои диэлектрические, возникает МП. С дру-

го коэффициента отражения при нормальном паде-

гой стороны, в длинноволновом пределе получены

нии. То есть при падении под углом параметры мо-

материальные уравнения для такой одноосной ди-

гут измениться. Результаты справедливы в предпо-

электрической ИС: εe⊥ = (ε1d1 + ε2d2) /d, ε-1e∥ =

(

)

ложении k₀a ≪ 1, где a — радиус шариков, тогда как

ε-11d₁ + ε-12d₂

/d, а также аналогичная форму-

для гомогенизации надо использовать более жесткое

ла для МП в случае ИС с магнитными слоями. Та-

условие k₀d ≪ 1, где a ≪ d, d — период. Если это

кие 1D-метаматериалы описываются диагональны-

жесткое условие выполнено, эффективная ДП мало

ми тензорами, в которые входят проницаемости сло-

отличается от проницаемостей основы, а МП умень-

ев и их толщины. Если магнитных свойств у сло-

шается на порядок и более. Таким образом, условие

ев нет, то μ1 = μ2 = 1 и μe = 1, т.е. нет их и

μ^′e ≫ 1 есть результат сильного взаимовлияния МА

у ИС. Однако такой метаматериал можно описать

(влияния внутреннего поля), при этом модель ста-

моделью, учитывающей магнитные свойства, введя

новится менее строгой из-за пространственной дис-

два скалярных параметра ε_e и μ_e вместо тензора ε_e.

персии и других факторов. Основное ограничение

Например, можно потребовать, чтобы коэффициент

модели при увеличении размеров — исчезновение

отражения R от такой среды был равен коэффици-

сферической симметрии, используемой при выводе

енту отражения от изотропной эффективной среды.

формул.

Коэффициент отражения зависит от угла падения,

В работе Рытова [19] рассмотрена плоскослои-

т. е. от вектора k в среде, поэтому подход позволя-

стая периодическая среда. Для нее найдены ска-

ет учесть ПД. Его можно распространить на струк-

лярные ДП и МП, описывающие волны двух попе-

туру типа мелкомасштабный период в крупномас-

речных направлений. Условия получены путем при-

штабном периоде [20]. Если структура имеет боль-

равнивания волнового сопротивления и коэффици-

шое число N периодов по два диэлектрических слоя

ента замедления своим эффективным значениям:

в периоде, то ее коэффициент отражения и коэффи-

197

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

циент прохождения рассчитываются точно. Удобно

Для исследования свойств метаматериалов при-

использовать матрицу передачи. Описывая струк-

менялись коммерческие пакеты программ, позволя-

туру как однородную с двумя параметрами ε_e и μ_e,

ющие рассчитывать конфигурационно сложные пе-

мы можем их связать с параметрами слоев. Такая

риодические микроструктуры. Такие пакеты могут

модель будет достаточно правильно описывать ди-

демонстрировать решения, иллюстрирующие отри-

фракцию на плоскослоистых структурах при дру-

цательную рефракцию (ОР) или преломление. Од-

гих больших значениях N, т. е. при других толщи-

нако ОР никак не означает, что ε < 0 и μ < 0, и

нах всей структуры. Однако она не применима для

даже ничего не говорит о том, что какие-то компо-

описания волн внутри. Более подходит модель на ос-

ненты тензоров отрицательны. Обычно ОР связы-

нове тензора ДП одноосного кристалла. Приведен-

вают с объемными ОВ. Однако ОВ и ОР — явления

ная выше модель, описывающая ε_e⊥ и ε-1e∥, является

разные и могут существовать независимо [1,4]. Час-

приближенной. Можно учесть ПД, сшив поля и вы-

то такие ИС называют левыми, метаматериалами

числив поляризацию под воздействием плоской вол-

с отрицательной групповой скоростью или с отри-

ны направления k [1]:

цательным показателем преломления. Отрицатель-

ным показатель преломления быть не может [21].

⌊

(

)

ε_xx = k-10

k2x + k21zd-2

d²¹ + d₁d₂

Это понятие введено в оптике на заре ее развития и

(

)⌋

применимо для изотропных сред, когда длина вол-

+ k22zd-2

d²² + d₁d₂

ны в 10⁵ раз и более превышает размеры молекул.

Для ИС это соотношение как минимум на три по-

рядка меньше, и для них ДП ε(ω, k) есть тензор,

εzz = k⁰²k^x ×

{

(

)

(

)}-1

причем зависящий от волнового вектора k [4, 6-8].

k21zd-2

d²¹+d1d2ε2/ε1

+k22zd-2

d22

+d1d2ε1/ε2

× 1-

Объемные волны в таком кристалле удовлетворя-

k^x+k21zd-2 (d²¹+d1d2) +k2zd-2 (d²+d1d2)

ют уравнению Френеля в общем случае четверто-

го порядка по k и шестого по волновому числу k₀.

Здесь обозначены ДП ε_n и толщина слоев d_n,

Максимальное количество таких волн — четыре, а

k2nz

= k²⁰ε_n - k2x, d

= d₁ + d₂. Интересно от-

вопрос о прямой или обратной объемной волне в

метить, что результат этих формул в ряде част-

недиссипативном фотонном кристалле решается по-

ных случаев совпадает с результатом приведен-

средством определения угла между векторами k и

ных выше. Нормированное к импедансу Z₀

v_g = ∇_kω (ε, k). Тип волны зависит от направления

√

μ₀/ε0 волновое сопротивление необыкновен-

групповой скорости (нормали к поверхности изо-

ной волны в таком 1D-фотонном кристалле равно

частот) по отношению к вектору k. В диссипатив-

√

Z₀ = E_x/H_y =

ном кристалле вместо v_g следует использовать век-

x - (k^x/k⁰)/(εxxεzz). В вакууме

√

тор Пойнтинга S = Re (E × H^∗)/2. В общем слу-

= Z₀

1 - k2x/k²⁰. Коэффициент отражения R

чае без бианизотропии поверхность изочастот опре-

= (Z/Z₀ - 1) / (Z/Z₀ + 1) связываем с тензором ДП,

деляется уравнением ω (ε(ω, k) , μ (ω, k) , k) = const,

что дает еще одну возможность теоретического

в котором учтена ПД. Для гиперболического мета-

и экспериментального определения его компонент.

материала это приводит к ограниченности поверх-

С (ру)ой стороны, из уравнения Флоке - Блоха

ности изочастот [1]. Обладающий магнитным от-

cos

= (a₁₁ + a₂₂)/2 можно определить блохов-

кликом бианизотропный метаматериал описывает-

ское волновое число

kи импеданс E-волны в направ-

ся еще двумя тензорами кросс-поляризации [14-16].

(

)

√

лении z: Z =

μ₀/ε₀k₀k/

k² + k²

. Здесь a_nn —

Он характеризуется сильной ПД. Однако он также

элементы матрицы передачи одного периода из двух

может быть описан только одним тензором ДП, но

слоев. Отличие двух импедансов в том, что первый

такой тензор в общем случае не может быть при-

построен на основе гомогенизации в виде уравнения

веден к диагональному виду. Для получения отри-

Френеля, а второй соответствует строгому решению

цательной рефракции важно, под каким углом вол-

уравнения Флоке - Блоха. Во вторую формулу для

на падает и какова оптическая плотность среды па-

коэффициента отражения эффективные проницае-

дения [4]. При падении из среды с ДП ε касатель-

мости не входят. Приравнивая оба коэффициента

ная к поверхности компонента вектора

k сохраня-

отражения, получим формулу для гомогенизации:

ется, а сам вектор удовлетворяет уравнению Френе-

ля

k₌ k²⁰ ε. Внутри кристалла вектор k удовлетво-

√

k2x/k²⁰

k₀

√ε_xxk

ряет более сложному уравнению Френеля. Изменяя

ε_zz

угол падения, мы изменяем положение конца век-

k² + k²

198

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

тора на поверхности изочастот. Ввести один пока-

ны быть симметричными, т. е. иметь одинаковые ди-

затель преломления для таких кристаллов нельзя,

польные моменты (коэффициенты поляризуемости)

даже если взять его тензорным. Да и вводить его

вдоль каждой из осей. Диссипация должна быть ма-

нецелесообразно: в электродинамике вводится век-

лой: |k^′′| / |k^′| ≪ 1, |n^′′/n^′| ≪ 1, где k = k^′ - ik^′′ —

тор k и волновые матричные импедансы, поскольку

волновой вектор, n =

√ε_eμ_e = n^′ - in^′′ — получен-

гомогенизированные (макроскопические) уравнения

ный гомогенизацией показатель преломления (заме-

Максвелла принимают вид

тим, что можно ввести вектор n^′ = k^′/k₀, модуль

которого определяет замедление в направлении k^′).

k × H = -Z₀k₀εe (ω,k)E,

ПД и биизотропия должны быть пренебрежимо ма-

k × E = -Z-10k₀μe (ω,k)H,

лыми. Последнее соответствует тому, что метама-

а импедансы определяются из них, когда определе-

териал следует рассматривать вдали от резонансов

на зависимость k(ω) для конкретной дисперсионной

и запрещенных зон при малых |k| и больших по

ветви и получены эффективные параметры ε_e (ω, k)

сравнению с размерами частиц длинах волн Λ. Пре-

μ_e (ω, k). В диссипативных ИС отрицательная

небрежение биизотропией означает, что вторичные

групповая скорость и объемные ОВ — понятия

электрические поля, создаваемые за счет поляриза-

не тождественные. При диссипации k-пространство

ции каждой частицы, не вносят вклад в ее магнит-

уже не является трехмерным: оно шестимерное ком-

ную поляризацию и в магнитную поляризацию ее

плексное, что не позволяет определить групповую

соседей. Это же относится и к магнитным полям.

скорость как градиент скаляра. Поэтому скорость

Еще раз подчеркнем, что использование в уравне-

движения энергии и групповая скорость — так-

ниях Максвелла скалярных величин ε < 0 и μ < 0

же понятия различные. В металлических фотонных

дает решение в виде плоской волны, для которой

кристаллах в области плазмонного резонанса |ε^′| ∼

k·S < 0, что не требует введения ни отрицательного

∼ ε^′′, а сделать ε^′′/ |ε^′| < 10-3 при комнатной тем-

показателя преломления, ни отрицательной группо-

пературе практически невозможно. Поверхностные

вой скорости. В монохроматической волне нет груп-

волны также могут быть обратными [22], и для это-

пы волн, поэтому нет и основания для ее введения.

го не требуется наличие магнитных свойств у сре-

ды, хотя есть и обратные магнитостатические вол-

ны [23]. В отличие от наведенного магнетизма в ИС

3. ПОЧТИ ИЗОТРОПНЫЕ

анизотропный магнитный отклик в ферритах про-

МЕТАМАТЕРИАЛЫ

является в СВЧ-диапазоне в присутствии внешнего

магнитного поля, что связано с ограниченностью на-

Электрический диполь длины l хорошо модели-

магниченности насыщения. Поверхностные ОВ мо-

руется малым тонким металлическим цилиндром

гут испытывать отрицательное преломление на ме-

малого радиуса r ≪ l. Изотропный диполь пред-

таповерхностях.

ставляет собой три таких взаимно перпендикуляр-

Рассмотренные выше понятия требуют коррект-

ных цилиндра (рис. 2a). Его излучение изотропно.

ного использования. Часто ряд из них лишен фи-

Магнитный диполь представляет собой проволоч-

зического смысла. Так, формальная подстановка в

ную рамку. Изотропный магнитный диполь предста-

уравнения Максвелла однородных и изотропных ве-

вим как три скрещенные рамки (рис. 2). Поместим

личин ε_e = ε^′ -iε^′′ и μ_e = μ^′ -iμ^′′ при ε^′′ > 0 и μ^′′ > 0

два (магнитный и электрический) изотропных ди-

дает для показателя преломления n_e = n^′ - in^′′ зна-

поля в каждый узел кубической решетки (рис. 2a).

чение n^′ < 0, если взять ε^′ < 0 и μ^′ < 0. Цель дан-

Для устранения взаимного влияния диполей сосед-

ной работы — рассмотреть, возможно ли получить

них узлов предположим, что R ≪ a. Считаем l < R.

изотропную среду с ε^′ < 0 и μ^′ < 0. Ниже показа-

Даже в этом случае имеет место ближнепольное вза-

но, что это невозможно, что в некотором роде пере-

имное влияние диполей в узле: возбуждение элект-

кликается с выводами работы [24]. Волны в таком

рического диполя приводит к излучению, которое

метаматериале в низкочастотном пределе подчиня-

возбуждает токи в рамках, а возбуждение рамок,

ются уравнению Френеля k² = k²⁰ε_eμ_e. Построен-

в свою очередь, возбуждает электрические диполи.

ный как кристалл, он должен обладать следующи-

Поэтому такой кристалл не совсем удовлетворяет

ми свойствами. Решетка должна быть кубической с

нашим требованиям. Он, скорее, будет биизотроп-

периодом a ≪ Λ, где Λ — минимальная внутренняя

ным. В нем каждый излучатель в узле принадле-

длина волны, а λ = Λ√εeμe — длина волны в ва-

жит восьми кубическим ячейкам. Поскольку в ку-

кууме. Включенные в узлы решетки частицы долж-

бе их восемь, одному кубу принадлежит один из-

199

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

Рис. 2. Фотонный кристалл с магнитными и электрически-

ми диполями в кубической решетке

лучатель. Его можно расположить в центре куба.

Тогда надо задать разложение токов в проволоч-

ках, вычислить для них поля, используя периодиче-

скую функцию Грина (ФГ) данного фотонного кри-

сталла и наложить граничные условия на проволоч-

ках. Это позволяет сформулировать дисперсионное

Рис. 3. Нормированная глубина проникновения в серебро

уравнение. Его будем называть микроскопическим.

в зависимости от обратной длины волны

Оно позволяет построить точную зонную структу-

ру и микроскопическую поверхность изочастот. Для

ее построения нужно дисперсионное уравнение, в

которое входят только поля, частота и волновой

приводить к полному проникновению поля в шари-

вектор k. Теперь можно выполнить гомогенизацию,

ки или проволочки. Для наноразмерных металли-

определить средние поля, индукции и поляризации,

ческих структур оно возможно в широких частот-

ных диапазонах, включая оптические и УФ-области.

что позволяет определить усредненные (эффектив-

ные) гомогенные в общем случае тензорные прони-

На рис. 3 приведена кривая нормированной к длине

цаемости. Им соответствует уравнение Френеля для

волны глубины проникновения в серебро в зависи-

фотонного кристалла, в которое входят величины

мости от обратной длины волны в метрах. Исполь-

зована формула Друде - Лоренца для объемного об-

εe, μe, k, k0. При фиксации k0 оно дает макроско-

пическую поверхность изочастот и может считаться

разца. Для наноразмерных пленок и квантовых ни-

макроскопическим дисперсионным уравнением.

тей эти результаты являются приближенными, по-

На рис. 2б представлена ИС, которая более со-

скольку вычисление требует квантовых подходов.

ответствует требованиям изотропности. Очевидно,

Различные резонансные излучатели типа разо-

три кольцевых рамки можно заменить металличес-

мкнутых колец, двойных разомкнутых колец и т.п.

ким шаром. Кольцевые токи на поверхности ша-

не являются изотропными. Фотонные кристаллы на

ра наводят магнитную поляризацию. Это низкочас-

их основе не являются изотропными средами. В об-

тотные токи. Шар также может представлять со-

щем случае метаматериалы с такими МА являют-

бой электрический диполь, квадруполь, октуполь и

ся бианизотропными: электрическое поле наводит в

т.п. [25]. Мультипольные резонансные токи соответ-

них как электрическую, так и магнитную поляри-

ствуют более высоким частотам, хотя они существу-

зации, а магнитное поле — соответственно магнит-

ют на любых частотах в силу разложений полного

ную и электрическую поляризации. Такие ИС могут

тока по мультиполям, т.е. по производным полино-

демонстрировать ОР, но никак не являются струк-

мов Лежандра. Однако ориентацию оси z можно вы-

турами с ε < 0 и μ < 0 или с n < 0. Их стро-

брать произвольно, поэтому это изотропные излу-

гий анализ весьма сложен. Другие почти изотроп-

чатели. Часто токи считают поверхностными. Для

ные структуры приведены на рис. 4. Они более про-

этого глубина проникновения должна быть суще-

стые для анализа и гомогенизации. Для металличе-

ственно меньше размеров МА. Говоря далее о токах,

ских шариков (рис. 4a) имеет место диамагнетизм

мы всегда подразумеваем их плотности, за исключе-

и малая эффективная ДП [26]. Металлические ци-

нием отдельно оговоренных случаев. Слишком низ-

линдры (рис. 4б) демонстрируют изотропную элек-

кие частоты или слишком малые размеры могут

трическую и магнитную поляризации. Диэлектри-

200

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

ву, что снижает эффект. Во-вторых, мода отдельно-

го резонатора и поле резонатора в фотонном кри-

сталле при движении волны в нем совершенно раз-

личны. Волна заданного направления k в кристал-

ле из таких МА соответствует определенной точке

на дисперсионной поверхности, и изменение часто-

ты ведет к изменению k. От этого зависит направ-

ление поляризации. В силу открытости сферическо-

го резонатора его моды являются квазисобственны-

ми с комплексными резонансными частотами. Они

не образуют полную систему функций, а МА об-

менивается энергией со всеми соседями, включая

удаленные. Поле внутри МА при возбуждении рас-

сматриваемой волной состоит из бесконечного на-

бора таких мод и функций непрерывного спектра,

т. е. нельзя выделить одну моду. В-третьих, рассмат-

риваемый кристалл обладает сильной ПД, т. е. не

является изотропным, поскольку на высоких резо-

Рис. 4. Элементарные ячейки кубических фотонных кри-

нансных частотах возбуждение шара зависит от на-

сталлов с шаровыми (а, г), цилиндрическими (б), шаровы-

правления k, а отклик не является локальным. В-

ми и кольцевыми (в) включениями МА в диэлектрическую

четвертых, в диэлектрических структурах магнит-

основу (матрицу)

ные свойства можно не вводить, описывая их только

тензором ДП. В-пятых, даже если из каких-то сооб-

ражений вводится магнитная поляризация, необхо-

ческие шары и металлические кольца (рис. 4в) де-

димо доказывать, что какие-то ее компоненты на-

монстрируют изотропную поляризацию обоих ти-

ходятся в противофазе к магнитному полю, чтобы

пов. Структура (рис. 4г) из двух сдвинутых пря-

утверждать об отрицательности компонент МП.

моугольных решеток с различными диэлектриче-

скими шарами в рассматриваемом низкочастотном

4. ИНТЕГРАЛЬНОЕ И ДИСПЕРСИОННОЕ

приближении может создавать только изотропную

УРАВНЕНИЯ

электрическую поляризацию. Такой метаматериал

представляет собой две сдвинутые на полпериода

Рассмотрим фотонный кристалл с кубическими

одинаковые кубические решетки с разными шаро-

выми МА в них. Замена металлических цилиндров

ячейками с размером граней a, заполненными ди-

на диэлектрические возможна. Замена шаров на

электрической средой с ДП ε, и диэлектрические

тела или МА в ней объема V с поверхностью S. Об-

симметрично расположенные кубики также возмож-

на. Кубики существенно удобнее для моделирова-

ласть может быть многосвязной, т. е. объем может

состоять из суммы объемов, ограниченных не кон-

ния. Но технологически кубики невозможно распо-

ложить в узлах решетки без случайных поворотов.

тактирующими поверхностями. Комплексную ДП

ε (ω) = ε^′ (ω) - iε^′ (ω) МА считаем не зависящей от

В ряде работ было предложено использовать ИС

рис. 4г как метаматериал с одновременно отрица-

координаты r внутри тела, а границы МА резкими,

т. е. ДП при переходе через границу терпит скачок:

тельными ε и μ [27-31]. Приводилось такое объясне-

ние. Пусть совпадают резонансные частоты каких-

ε^- (ω, r) = ε⁺ (ω, r) = ε, r ∈ S. При этом на по-

верхности S терпит скачок и нормальная компонен-

либо E-мод и H-мод для двух различных шари-

ков. Тогда они дают вклад в поляризацию как элек-

та электрического поля: E^-n ε = E+n ε, а на грани-

трическую, так и магнитную. Чуть сдвигая частоту

це имеются источники — наведенная поверхностная

плотность связанных зарядов. Очевидно, что ее ин-

вверх, мы как бы должны получить поляризации,

сдвинутые по фазе относительно полей. Здесь мож-

теграл по всей поверхности, т. е. полный связанный

заряд, равен нулю. В такой постановке задачу удоб-

но возразить так. Во-первых, все резонансы высо-

кочастотные. Известные изотропные диэлектрики с

но решать методом объемных интегральных урав-

нений с использованием периодической скалярной

малыми потерями не обладают очень большими ДП.

МА необходимо включить в диэлектрическую осно-

ФГ [1]:

201

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

горитмизации. Мы будем использовать объемное ин-

^G(r) =

a³

тегральное уравнение

(

⌊

⌋)

∑

exp

-i

kxkx +^kyly +^kzmz

E (r) (2 + ε (r) /ε) /3 =

(1)

+k2yl +k2zm - k²⁰ε

∫

k,l,m=-∞

= p.v.

L Ĝ (r - r^′)(ε(r^′)/ε- 1)E(r^′) d³r^′.

(2)

Здесь

kxk = kx + 2kπ/a,

kyl = ky + 2lπ/a,

kzm = kz +

+ 2mπ/a, что и определяет периодичность по k^′ с

точностью до пространственных гармоник. В слу-

Здесь обозначен оператор

L ≡ k

I + ∇ ⊗ ∇. Ин-

чае учета потерь k = k^′ - ik^′′. ФГ (1) удовлетворяет

теграл в (2) берется как главное значение по Ко-

уравнению

ши. Зависимость ФГ от k и k опущена, но как раз

она определяет дисперсионное уравнение. В пра-

(

)

∇² + k²⁰ε

^G(r - r^′) = -δ (r - r^′)exp(-ik · r^′) .

вой части (2) не зависящий внутри объема одно-

родных частиц множитель κ = ε/ε - 1 выносим

Диэлектрическое тело c ДП ε в основе c ДП ε

из-под знака интеграла. В силу ФГ (1) это интег-

создает плотность тока поляризации J_p (r)

ральное уравнение достаточно решить только в од-

= iωε₀ (ε - ε)E (r). Ток поляризации является до-

ной (нулевой) ячейке периодичности и только внут-

полнительным к току смещения в основе. Отсю-

ри частиц, где ε (r) = ε. Поэтому множитель у по-

да нетрудно видеть, что в среде с отрицатель-

ля в левой части есть 1 + κ/3. ФГ обеспечивает

ной ДП имеет место сдвиг по фазе на

-π/2

правильное взаимовлияние всего бесконечного ан-

между током поляризации и полем, тогда как в

самбля МА. В (2) можно использовать соотноше-

обычном диэлектрике этот сдвиг равен π/2. Та-

ние

L G (r - r^′)

= L^′ G (r - r^′), где штрих означа-

кой подход применим для всех частот и любых

ет дифференцирование по штрихованным координа-

рассматриваемых МА. Для сверхнизких частот

там (точке истока). Поле E будем аппроксимировать

(

)

J_p (r)

ε₀ω2p/ω_c

E (r)

= σ₀E (r), где обозна-

объемными кубическими кусочно-постоянными ко-

чена проводимость на постоянном токе. На низ-

нечными элементами, нумеруя их одномерным ин-

ких частотах это проводимость Друде. В оптике

дексом m. Пусть таких элементов M. Метод Га-

существенен член Лоренца. Для тонких и длин-

леркина для

(2) приводит к однородной систе-

ных проволочных структур важна только продоль-

ме линейных алгебраических уравнений (СЛАУ)

ная компонента тока, что упрощает расчет. Инте-

(1 + κ/3)IE = κĞE размерности 3M, где

I — еди-

грал от ФГ (1) с плотностью J_p (r) дает вектор-

ничная матрица размерности 3M, E — вектор-стол-

потенциал A (r), дифференцируя который прихо-

бец размерности 3M, состоящий из вектор-столбцов

дим к объемному интегральному уравнению в ви-

E^x, E^y и E^z размерности M. Значения E^xm, E^ym и E^zm

(

)

де E (r) = (iωε₀ ε)-1

∇∇ · A (r) + k²⁰ εA (r)

. Далее

соответствуют значениям соответствующих компо-

√

обозначим k = k₀

ε. Кроме этого имеем магнитное

нент электрического поля в элементе с номером m.

поле H (r) = ∇ × A (r). Можно записать несколько

Матрица

Ğ размерности 3M имеет блочную струк-

Ğ^αβ

уравнений для поля E, для поля H, для их ком-

туру

, где α, β принимают значения x, y, z.

mm^′

бинации, а также уравнения, нагруженные поверх-

Дисперсионное уравнение имеет вид равенства ну-

ностными интегралами по границе S, и без них [1].

лю определителя указанной СЛАУ. Находить корни

В общем случае они имеют форму интегродиффе-

такого уравнения весьма сложно. Удобнее исполь-

ренциальных уравнений, поскольку искомые вели-

зовать стационарный квадратичный функционал,

чины находятся как под знаком интеграла, так и под

получающийся умножением (2) скалярно на E^∗(r)

знаком производных. Использование поверхностных

и интегрированием по объему. Его удобно разре-

интегралов, как видно из рис. 2, 4, неудобно для ал-

шить относительно k²:

∫

(1 + k/3)

|E (r)|² d³r - κ · p.v.

E (r) ∇ ⊗ ∇G (r - r^′) d³r^′ d³r

V V

k² =

∫ ∫

(3)

E^∗ G (r - r^′)E (r^′) d³r^′ d³r

V V

202

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

Правая часть (3) нелинейно зависит от k, поэтому

из условий экстремума функционала. Часто, напри-

уравнение (3) следует решать совместно с уравне-

мер для тонких и длинных цилиндров, можно пре-

нием (2) итерационно. Удобно сначала найти при-

небречь какими-то компонентами поля (радиальной

ближенное решение дисперсионного уравнения, за-

и/или азимутальной).

дать одну из амплитуд, например, Ex1, что задает

амплитуду волны, и решить СЛАУ прямыми мето-

дами. Это позволяет уточнить решение дисперсион-

5. ГОМОГЕНИЗАЦИЯ

ного уравнения и далее использовать итерации. Ес-

ли по каждой из координат разбить частицу на 10

Гомогенизация — процедура (в общем случае

конечных элементов, то всего получим 10³ элемен-

неоднозначная) получения эффективных парамет-

тов, а размерность задачи будет 3·10³. Для матрицы

ров однородных (гомогенных) сред, эквивалент-

СЛАУ и функционала (3) можно написать квадра-

ных электродинамически в некотором смысле исход-

турные формулы. Диагональные элементы при вы-

ным ИС, основанная на решении обратных задач. В

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

общем случае для нее нужны некоторые процеду-

можно приравнять нулю или использовать другие

ры усреднения. Зависимость методов усреднения —

приближенные оценки.

одна из причин неоднозначности. Другая причина

В принципе дисперсионное уравнение позволя-

заключается в неоднозначности модели гомогенной

ет при задании k определить ω. Без учета дис-

среды. Третья причина заключается в разных мето-

сипации k можно задавать произвольно. Фикси-

дах гомогенизации. Один из первых методов, при-

руя k₀, можно построить поверхность изочастот в

мененный в теории искусственных диэлектриков, —

k-пространстве. Величина v_g = ∇ω (k) есть группо-

использование параметров дифракции [17]. Для вы-

вая скорость, показывающая направление движения

полнения условий теоремы погашения ИС должна

энергии. В случае диссипации групповую скорость

иметь много периодов по каждому измерению. Дан-

использовать нельзя, следует вычислять усреднен-

ный подход электродинамически сложный, особен-

ный по ячейке вектор Пойнтинга 〈S〉. Он может

но для 3D-задач. Второй метод основан на сравне-

не совпадать по направлению с v_g. При этом век-

ниях расчетов дисперсии по электродинамическим

тор k = k^′ - ik^′′ является комплексным. Подлежа-

моделям (например, по уравнениям типа (2), (3)) с

щих определению скалярных компонент становится

дисперсией на основе уравнения Френеля. Он отно-

больше (шесть) и одного комплексного дисперсион-

сительно простой, если получено микроскопическое

ного уравнения уже недостаточно. k-пространство

дисперсионное уравнение. Третий способ, который

также становится шестимерным. Необходимо вы-

будем использовать здесь, это гомогенизация на ос-

числять вектор Пойнтинга S = Re (E × H^∗) /2 и ис-

нове вычисления усредненных по ячейке периодич-

пользовать условие S/ |S| = k^′′/ |k^′′|. Его смысл в

ности поляризацией для волны, определяемой вол-

том, что направление затухания волны k^′′ совпада-

новым числом k₀ и волновым вектором k. Концеп-

ет с направлением движения энергии.

туально он соответствует определению параметров

макроскопической электродинамики путем усредне-

Поскольку дисперсионное уравнение в виде ра-

венства нулю определителя большого порядка очень

ния по физически бесконечно малому объему. Од-

нако размер ячейки периодичности a не обязатель-

сложное, более простой приближенный подход мо-

жет быть основан на использовании уравнения в ви-

но мал по сравнению с длиной волны. Во всяком

де функционала, построенного на основе (2). Для

случае, при строгом решении уравнений Максвел-

этого умножим (2) на E^∗ и проинтегрируем по объ-

ла, усреднения возможны при любых соотношени-

ему или МА. Приближенное уравнение получается,

ях, если только не накладывается условие a ≪ λ.

если выбрать электрическое поле (или ток) из физи-

Это условие важно для введения локальных про-

ческих соображений. Для проволочных МА, когда

ницаемостей. Введем усредненные по ячейке поля,

глубина проникновения меньше их радиуса, вмес-

поляризации и индукции, например,

то поля удобно использовать осевой ток, а гранич-

∫

ные условия налагать на поверхности проволочек.

〈E_x〉 =

E_x (r) d³r.

a³

Это приводит к простым дисперсионным уравнени-

ям [26]. Для МА в виде однородного диэлектриче-

ского цилиндра или шара удобно использовать яв-

Электрическую поляризацию запишем как

ные низкочастотные решения уравнения Гельмголь-

(

)

ца внутри и определять неизвестные коэффициенты

〈P^e〉 = ε₀

εe -

〈E〉 + c-1

ξ 〈H〉 .

203

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

Здесь мы рассмотрели бианизотропную модель.

писать соотношения только для одной компоненты.

Аналогично для магнитной поляризации

Имеем

(

)

〈P^m〉 = μ₀

μ_e -

〈H〉 + c-1

ζ 〈E〉 .

ε_e = ε+ ε-10 〈P^ex〉/ 〈E_x〉 ,

(4)

(

)

μ_e = 1 + μ-10 〈P^mx〉/ 〈H_x〉.

Введем намагниченность M =

μ_e -

I H. Запи-

сывая другие усредненные компоненты, получаем

Получив решение задачи для E, поле h (r)

СЛАУ для определения эффективных параметров.

= ∇ × A(r) находится как H =

AE, где введен-

Подход годится для анизотропной и бианизотроп-

ная матрица определяется интегродифференциаль-

ной моделей сред с учетом соотношений симметрии

ным оператором получения магнитного поля. Рас-

для материальных параметров. При большом чис-

смотрим вычисление поляризаций:

ле неизвестных можно написать переопределенную

∫

СЛАУ и решать ее методом регуляризации. В случае

〈P^e〉 =

ρ (r) r d³r,

изотропной модели

ξ =

ζ= 0. В случае биизотроп-

a³

ной модели эти параметры — ненулевые псевдоска-

(5)

ляры. В случае диэлектрической анизотропной сре-

∫

ды имеем три уравнения:

〈M〉 =

r × Jp (r) d³r.

a³

〈P^ex〉 = ε₀ ((ε_exx - ε) 〈E_x〉 + ε_exy 〈E_y〉 + ε_exz 〈E_z〉)

Здесь имеется плотность заряда ρ (r) = i∇ · J_p/ω =

и два аналогичных для других компонент. При пре-

= -ε₀∇ · [(ε - ε)E]. Поскольку мы рассматри-

небрежении диссипацией имеют место условия сим-

ваем МА из однородного диэлектрика, внутри

метрии Онзагера - Казимира ε_eαβ = ε_eβα, т. е. тен-

него ∇ · E. Поэтому следует вычислить поверхност-

зор является симметричным. При диссипации тен-

ную дивергенцию, связанную с наведенной поверх-

зор эрмитовым не будет, но будет выполнено условие

ε′eαβ = ε′eβα. Из трех уравнений определяются три

ностной плотностью заряда ρ_S

= ε₀E_ν (1 - ε/ε).

Здесь E_ν — внешняя нормальная координата. По-

величины, если учесть условия симметрии. Уравне-

ния упрощаются при приведении тензора ДП к глав-

скольку ∇ · (εE) = 0, имеем ε∇ · E = - (ε/ε) E · ∇ε,

∇ε = (ε - ε), и плотность заряда пропорциональна

ным осям. Итак, по трем поляризациям можно опре-

дельта-функции от нормальной координаты. Поэто-

делить три компоненты тензора ДП. Если есть толь-

му первый интеграл в (5) становится поверхност-

ко два тензора ε_e и μ_e и их можно одновременно при-

вести к диагональному виду, то шесть в общем слу-

ным:

∮

чае комплексных компонент можно определить из

〈P^e〉 =

ρ_S (r)rd²r.

a³

шести независимых комплексных скалярных урав-

нений Максвелла. Однако возможность такого при-

Обозначим χ = ε/ε- 1. Уравнение (2) можно запи-

ведения зависит от структуры ИС. Например, пусть

сать в следующей форме: E =

^Ĝ (χE),

Ĝ =

LG(),

у нас имеется ИС с магнитными и электрическими

диполями в узлах двух одинаковых вложенных ку-

где

G — интегральный оператор с ядром

^G,

L =

бических решеток. Мы можем изготовить ИС, сдви-

= ∇ ⊗ ∇ + k² — дифференциальный оператор. По-

(

)

нув решетки друг относительно друга на произволь-

скольку ∇ · ∇ ⊗ ∇ = ∇ · ∇² и

∇² + k²

G₌ -δ (r),

ное расстояние. Формально мы также можем повер-

нетрудно видеть, что действие оператора диверген-

нуть одну решетку относительно другой на три про-

ции слева на интегральное уравнение дает тождест-

извольных угла относительно разных осей. Навряд

во. Решая это уравнение, следует определять нор-

ли следует реализовывать такой мысленный экспе-

мальную компоненту электрического поля и поверх-

римент, но это отличает ИС от природных сред,

ностную плотность заряда.

где симметрия существенна. В случае только ди-

Рассмотрим ряд простых случаев. Пусть имеем

электрической модели ИС без диссипации условия

фотонный кристалл из проволочных колец, ориен-

Онзагера - Казимира выполняются и в общем слу-

тированных вдоль оси z. Для него ε_e = 1, μ_ezz = μ,

чае тензор приводится к диагональному виду. Мо-

μ_exx = μ_eyy = 1. Имеем одну компоненту электри-

жет оказаться, что более удобна система координат,

ческого тока J_ϕ = σ (ω) E_ϕ. Пусть для простоты ра-

в которой достаточно определить шесть его компо-

диус проволочки меньше толщины скин-слоя, а ток

нент из шести уравнений Максвелла, а не три. Для

распределен равномерно. Соленоидальный ток со-

интересующей нас изотропной ИС достаточно за-

здает две компоненты вектор-потенциала [25]:

204

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

∫

а для ФГ запишем выражения

A_ρ = J_ϕ sin(ϕ - ϕ^′)

^G(r, r^′) d³r^′,

∫

G_xlm = exp(-ilk_xa - imk_za)G_lm,

A_ϕ = J_ϕ cos(ϕ - ϕ^′)

^G(r, r^′) d³r^′.

G_ylm = exp(-ilk_ya - imk_za)G_lm.

ФГ G_xlm соответствует источник, расположенный в

Интегрирование проводится по объему кольца, а ток

точке x = la, y = 0, z = ma, а ФГ G_ylm — источ-

можно вынести из-под интеграла, поскольку он по-

ник, расположенный в точке x = 0, y = la, z = ma.

стоянен. Считаем проволочки тонкими, а поле пол-

Однако поля от всех источников рассматриваются в

ностью приникающим в них. Нам удобнее использо-

нулевой ячейке. Поэтому с помощью суммирования

вать обычную ФГ

∑

G₌

exp(-ia [mk_z + l (k_x + k_y)])G_lm

G₀ (r, r^′) = |4π (r - r^′)|-1 exp(-ik₀ |r - r^′|),

m=-∞ l=-∞

получаем периодическую ФГ, действующую в нуле-

а затем ее следует периодически продолжить. Име-

вой ячейке, связанной с началом координат. Обозна-

ем

чим ψ_lm (k_x, k_y, k_z ) = a ⌊mk_z + l (k_x + k_y)⌋. Для ФГ

∇ · A = ρ-1∂_ρ(ρA_ρ) + ρ-1∂_ϕA_ϕ,

имеем [25]

(

√

)

∫

∞

∑

exp

κ² - k²⁰ |z - z^′ - ma| J_n (κρ)J_n (κ (ρ^′ + la))

G_lm =

exp(-in (φ - φ^′))

√

κ dκ.

(6)

4π

κ² - k²

n=-∞

Очевидно, что G₀ (r, r^′) = G₀₀ (r, r^′). Вычисляя ин-

∫

тегральное уравнение в нулевой ячейке, получаем

A_ϕ = 2

g1lm (ρ, z | ρ^′, z^′)ρ^′ dρ^′dz^′.

решение задачи. Чтобы не переходить в тороидаль-

ную систему координат, мы, оставаясь в цилиндри-

Мы воспользовались равенством J-1 (x) = -J₁ (x)

ческой, заменим круглое сечение проволоки радиу-

для функций Бесселя. Вектор-потенциал имеет од-

са r₀ на прямоугольное площади 4r²⁰, почти не изме-

ну компоненту, не зависящую от ϕ. Электрическое

няя результат. Считая поле постоянным, имеем

поле имеет три компоненты E_ρ = ∂_ρA_ϕ/ (iωε₀), E_ϕ =

∫

= k²⁰A_ϕ/ (iωε₀), E_z = ∂_zA_ϕ/ (iωε₀). Магнитное по-

sin(ϕ-ϕ^′)exp(-in (ϕ-ϕ^′)) dϕ^′ = i (δn,-1-δn,1),

ле имеет две компоненты H_ρ

= -∂_zA_ϕ и H_z

= ρ-1∂ρ (ρA_ϕ). В декартовой системе имеются все

шесть компонент. Компонентами E_ρ и E_z внутри

∫

проволочки можно пренебречь по сравнению с по-

cos(ϕ - ϕ^′)exp(-in (ϕ - ϕ^′)) dϕ^′ = δn,-1 + δn,1.

стоянной азимутальной компонентой поля. Внут-

ри проволочки электрическое поле E_ϕ равно сумме

всех продуцированных периодической ФГ откликов:

Обозначая интеграл в (6) как g_nim (ρ, z | ρ^′, z^′), по-

лучаем

E_ϕ =

∫

A_ρ = i (g-1lm (ρ, z | ρ^′, z^′) -

∑

exp(-ia [mk_z+l (k_x+k_y)]) E_ϕlm,

(7)

- g1lm (ρ,z | ρ^′,z^′))ρ′ dρ^′dz^′

= 0,

m=-∞ l=-∞

(

)

√

∫

∞

exp

κ²-k²⁰ |z-z^′-ma| J₁ (κρ)J₁ (κ (ρ^′+la))

πk²σ (ω)

E_ϕlm =

ρ^′dρ^′ dz^′

√

κ dκ.

(8)

iωε₀

κ² - k²

205

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

В (7) мы проинтегрировали по углу, поэтому ин-

⁽πk²⁰ |σ (ω)|)2 ×

тегрирование в (8) теперь выполняется по сечению

ωε₀

провода. Считая поле постоянным, умножим (7) на



ρ, проинтегрируем еще раз по сечению и сократим



∑



exp(-iψ_ml (k_x, k_y, k_z)) ×

на поле. Слева в (7) останется 4Rr²⁰. Правую часть



m=-∞ l=-∞

(8) рассмотрим отдельно для нулевого и ненулево-

√

го значений индекса m, обозначив α =

κ² - k²⁰.

× (2 - δm0) (2 - δl0) ×

В первом случае область интегрирования по штри-

хованной координате разбиваем на две. Результат



∫

∞



имеет вид

f_m (κ, k₀, r₀)I_l (κ, k₀, r₀)J₁ (κR)



√

κ dκ

(9)

r₀

κ² - k²⁰



f₀ (κ, k₀, r₀) =

Оно связывает частоту и волновой вектор. Этот ре-

∫r0

)

( √

зультат показывает, что в некоординатных задачах

exp

- κ²-k20 |z-z^′| dz dz^′ =

решения принимают сложный вид даже при силь-

−r0 -r0

(

)

ных упрощениях. В низкочастотном пределе α =

4r₀

1 - exp(-2αr₀)

= κ, и из (9) можно явно выразить k²⁰. Магнит-

2αr₀

ный момент для задачи принимает вид 〈M_z〉

= Jϕ (πr₀R)2 /a³. Для решения вопроса о знаке МП

Если величина α мала, то κ² ≈ k²⁰ и f = 4r²⁰. Во

необходимо вычислить усредненную компоненту по-

втором случае интегралы с положительными и рав-

ля:

ными по модулю отрицательными индексами равны

∫

〈H_z〉 =

dzH_z (r, k, k₀) .

f_m (κ, k₀, r₀) = f_-|m| (κ, k₀, r₀) =

a³

−a/2

-a/2

)₂

⁽sh(αr₀)

= -4 exp(-α|m| a)

Такие вычисления удобнее проводить в декартовой

системе координат. Знак МП определяет формула

Переходя к суммированию по положительным ин-

μ = 1 + 〈M_z〉/〈H_z〉. Если диссипации нет, намагни-

дексам, получаем 2fm (κ, k0, r0). Вычисляя инте-

ченность должна быть в противофазе к полю, т. е.

гралы по ρ^′, будем интегрировать по областям

отрицательной, и по модулю превышать единицу.

la - R - r₀ < ρ^′ < la - R + r₀, la + R - r₀ < ρ^′ < la +

Покажем, как эта задача решается в декартовой

+ R + r₀. Получаем

системе. Далее вместо (9) выведем новое диспер-

сионное уравнение, считая, рамку прямоугольной.

Кольцо радиуса R заменим прямоугольной рамкой

∫

с плечом b = R√π и сечением 4δ2 = πr20. Вынесем

Il (κ, r0, R) =

J₁ (κρ)ρ dρ +

из-под интеграла плотность J, направленную по ча-

la-R-r0

совой стенке вдоль осей x и y. Получаем компоненты

∫

вектор-потенциала и поля:

J₁ (κρ)ρ dρ ≈ 2r₀ [(la-R)J₁ (κ (la-R)) +

la+R-r0

-16iJ

A_x (r) =

+ +(la + R)J1 (κ(la + r))].

a³

(

⌊

⌋)

∑

exp

-i

kxkx +^kyly +^kzmz

Здесь R — радиус кольца. Мы воспользовались ма-

лостью радиуса проволочки и теоремой о среднем

+k2yl +k2zm - k²

k,l,m=-∞

значении интеграла. Имеем I_-l = I_l, поэтому сумми-

(

)

(

)

(

)

руем по положительным индексам, удваивая I_l. По-

sin

kzmδ

sin

kxkb/2 sin

kylδ

лучим аналогичный результат для комплексно-со-

пряженного дисперсионного уравнения. Учитывая,

kzm

kxk

kyl

(

))

что σ (ω) = iωε₀

ε_L - ω2P /

ω² - iωω_c

и перемно-

(

)

жая оба результата, получаем уравнение

× sin

kylb/2 ,

(10)

206

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

16iJ

Теперь надо усреднить поле (12). Усредняем следу-

A_y (r) =

a³

ющим образом:

(

⌊

⌋)

∑

exp

-i

kxkx +^kyly +

kzmz

+k2yl +k2zm - k²

(

)

∫

(

)

k,l,m=-∞

(

)

(

)

(

)

kxk

exp

-i^˜xkx dx =

a³

sin

kzmδ

sin

kxkb/2 sin

kylδ

−a/2

(

)

kzm

kxk

kyl

2 sin

k_xka/2

(

)

(

)

= sinc

kxka/2

× sin

kxkb/2 ,

(11)

kxka

Поэтому

4Jb²δ

H_z = ∂_xA_y (r) - ∂_yA_x (r) = -

a³

(

⌊

⌋)

4Jb²δ²

∑

exp

-i

kxkx +^kyly +

kzmz

〈H_z〉 = -

a³

(

)

(

)

(

)

+k2xk +k2xk - k²

k,l,m=-∞

∑

sinc

kxka/2 sinc

kyla/2 sinc

kzma/2

× F_mlk (k),

(12)

+k2yl +k2zm - k

k,l,m=-∞

(

)

(

)

× F_mlk (k).

(14)

F_mlk (k) = sinc

k_zmδ sinc

kxkb/2 ×

(

)[(

)₂

(

)

(

)₂

Функции sinc = sin (x) /x четные и быстро убыва-

× sinc

kylb/2

kyl

sinc

kylδ +

kxk

ют. В (14) существенный вклад вносит нулевой член

]

(

)

суммы, но строгий результат требует учета большо-

× sinc

kxkδ

(13)

го числа членов. Теперь запишем

χ_ezz = μ - 1 = -

(

)

(

)

(

)

(15)

∑

sinc

kxka/2 sinc

kyla/2 sinc

kzma/2

F_mlk (k)

+k2yl +k2zm - k²

k,l,m=-∞

В (15) следует подставлять решение микроскопичес-

при брэгговском резонансе в запрещенной зоне. Для

кого дисперсионного уравнения. Если же подста-

медленной волны максимальное замедление имеет

вить квадрат волнового числа из макроскопическо-

место при k_x = k_y = π/a, откуда получаем мак-

(

)

го уравнения Френеля k20 =

k2x + k2y

/μ_ezz + k2z,

симальное значение μ = 2π²/ (k₀a)² > 1. Частоту

то, задавая k, получаем неявное уравнение. Его ре-

максимального замедления можно найти из диспер-

шение симметрично по волновому вектору, μ (k) =

сионного уравнения. Выше частоты резонанса в за-

= μ(-k). Упрощение может быть достигнуто вы-

прещенной зоне μ < 0. Однако это так в предполо-

бором направления распространения волны. Так,

жении, что уравнение Френеля еще справедливо. Та-

выбирая k_x = k_y, k_z = 0, из уравнения (9) най-

кой подход для кристалла из металлических рамок

дем частоту и подставим в (15). Имеем F₀₀₀ (k_⊥) =

использован в работе [26]. Здесь имеет место анало-

= 2sinc(k_xδ) [k_xsinc(k_xb/2)]² и следующее нулевое

гия с кристаллооптикой, для которой сильная ПД и

приближение:

брэгговские резонансы наступают в рентгеновском

диапазоне. Используя же приближенную формулу с

2k2x - k²⁰

μ₍₀₎ = 1 -

учетом только одного члена ряда, получаем

2 [k_xsinc (k_xa/2) sinc (k_xb/2)]²

2π² - k₀a

Оно весьма грубое, что определяется медленной схо-

μ₍₀₎ = 1 -

≈1-

8sinc² (πb/ (2a))

димостью ряда (15). Для медленной волны знак вос-

приимчивости χ^ezz отрицательный. Из уравнения

Мы учли, что b ≪ a. В общем случае обозначим

Френеля следует 2k2x = μk20, поэтому в низкоча-

сумму в знаменателе (15) как

стотном пределе имеет место решение μ₍₀₎ = 1. Ес-

ли μ < 0, то волновые числа мнимые. Это возможно

S = S0 (k) +

S (k) +

S (-k) ,

207

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

S = S00m (k) + Sk00 (k) + S0l0 (k) + S0lm (k) +

Заменим металлическое кольцо тонким диэлект-

рическим. Пренебрежем поперечными полями по

+ Sk0m (k) + Skl0 (k) + Sklm (k),

сравнению с тороидальным. По кольцу идет ток по-

ляризации, который возбуждает такое же поле, как

где

и ток проводимости. Для металлического кольца нет

разницы: высокочастотный ток в той или иной сте-

sinc (k_xa/2) sinc (k_ya/2) sinc (k_z a/2)

S₀ (k, k₀) =

пени тоже есть ток поляризации. Для диэлектри-

k2x + k2y + k2z - k²

ческого кольца отсутствует ток на нулевой частоте.

В приведенных выше формулах для проводимости

×F₀₀₀ (k) ,

достаточно положить ω2P = 0. Весьма заманчиво по-

лучить магнитные свойства ИС на структурах ти-

S00m (k, k₀) =

па диэлектрические кольца, диэлектрические спи-

(

)

(

)

(

)

рали и вообще на МА с соленоидальными токами

∑

sinc

kxka/2 sinc

kyla/2 sinc

kzma/2

поляризации. Однако токи должны быть строго со-

k2x + k2y +k2zm - k²⁰

леноидальными. Для диэлектрического кольца ра-

m=1

диуса R должно выполняться условие k₀

√εR ≪ 1,

×F00m (k) ,

что для больших ДП приводит к весьма малым

радиусам и малым эффективным проницаемостям.

S0lm (k, k₀) =

ДП диэлектрического кольца, расположенного в ди-

(

)

(

)

электрической основе, должна сильно отличаться

∑

sinc (k_xa/2) sinc

kyla/2 sinc

kzma/2

от ее ДП. Диэлектрическое кольцо должно быть

k2x +k2yl +k2zm - k²⁰

очень тонким, чтобы можно было пренебречь по-

l,m=1

перечными полями и резонансами. Это приводит к

×F0lm (k) ,

малым эффективным параметрам. В отличие от ме-

S_klm (k, k₀) =

талла, на низких частотах ток поляризации мал, а

(

)

(

)

на высоких частотах приближение работает плохо.

∑

sinc (k_xa/2) sinc

kyla/2 sinc

kzma/2

Тем не менее, чисто формально можно использо-

k2x +k2yl +k2zm - k²⁰

вать полученные результаты для диэлектрических

k,l,m=1

колец. В вакууме такое кольцо является резонато-

×F_klm (k).

ром с H01δ-топом колебания (см. [29]). В фотон-

ном кристалле такой резонатор является источни-

Все остальные суммы вычисляются аналогично:

ком волны и в то же время возбуждается ей. По-

тройные, двойные и однократные, взятые по поло-

ляризация и ток кольца должны быть в противо-

жительным индексам. Тогда

фазе компоненте H_z. При отрицательной компонен-

те МП фотонный кристалл становится гиперболи-

μ(k, k₀) = 1 -

ческим метаматериалом. Рассматривая квадратную

S₀ (k) +

S (k) +

S (-k)

диэлектрическую рамку, можно использовать ком-

Имеем

S (k) =

S (-k), поэтому μ (-k, k₀) = μ (k, k₀).

поненты (10), (11) для получения дисперсионного

При k = 0 решением дисперсионного уравнения яв-

уравнения. Электрическое поле определяется через

(

) ˜G((ε-1)E),

ляется k₀ = 0, поэтому μ (0, 0) = 1. При малых раз-

вектор-потенциал как E =

k²⁰ + ∇∇·

мерах δ и b имеем μ (k, k₀) ≈ 1. Сильная ПД означа-

где A = iωε₀G ((ε - 1) E), а

G — интегральный опе-

ет, что нельзя получить зависимость μ (ω): резуль-

ратор с ядром (1). Использовать алгоритм с яд-

тат существенно зависит от k. Полученное матери-

ром ∇⊗∇G проблематично в силу неинтегрируемой

альное уравнение нелокально: оно относится ко всей

особенности. Считая поле соленоидальным внутри

области ячейки, а локальность наступает при λ ≫ a,

кольца и перенося оператор ∇· на (ε - 1) E, по-

т. е. при k² ≈ k₀, когда и μ (k, k₀) ≈ 1. Путем уве-

лучим поверхностные интегралы от поверхностной

личения размеров рамки можно получить значения

плотности заряда. Поскольку мы пренебрегли на

МП, отличные от единицы и даже отрицательные,

границе кольца нормальными компонентами поля,

но не факт, что реально достигается отрицательная

пренебрегли и поверхностными зарядами, и хоро-

МП в обычном понимании: в этом случае модель

шим приближением является интегральное уравне-

становится весьма приближенной, а результат зави-

ние E = k²⁰G ((ε - 1) E) = k20 (ε - 1)G (E). Диспер-

сит от k.

208

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

сионное уравнение на его основе запишем в виде

размеров указанная симметрия полей, вообще го-

воря, не выполняется. Она означает наличие двух

∫

электрических стенок, перпендикулярных плоско-

|E|² d³r = k²⁰ (ε - 1) ×

сти рамки, и может быть реализована при сдвиге

∫

фаз π на ячейку по обеим координатам, что может

× E∗ (r)

^G (r, r^′) E (r^′) d³rd³r^′.

(16)

выполняться в запрещенной зоне, которой соответ-

ствует полоса заграждения фильтра из конечного

В декартовой системе удобнее рассмотреть прямо-

образца. При дифракции на метаматериале в низ-

угольную рамку. Считая теперь поле постоянным и

кочастотном пределе реализуется сдвиг фаз, близ-

имеющим одну постоянную компоненту E_x или E_y

кий к нулю. Использование постоянных по значени-

в каждом плече, E_x (x, -b/2, z) = -E_x (x, b/2, z) =

ям компонент электрического поля внутри рамки —

= E_y (b/2,y,z) = -E_y (-b/2,y,z) = E, получаем

наиболее жесткое ограничение модели. Но именно

уравнение

оно позволило получить простой явный результат.

∑

4k²⁰δ²b (ε - 1)

|G_klm (k)|²

6. ОТРИЦАТЕЛЬНЫЕ ε^′e И μ^′e В МОДЕЛИ

a³

+k2yl +k2zm - k²

k,l,m=-∞

ЛОРЕНЦА

где

Классическая модель Лоренца для разреженного

(

)

(

)

газа осцилляторов (электрических диполей или ос-

|G_klm (k)|² = sinc²

kzmδ sinc²

kxkb/2 ×

цилляторов) с частотами ω0el дает выражение для

(

)[

(

)

(

)]

ДП [13]:

× sinc²

kylb/2

sin²

kxkb/2

+ sin²

kylb/2

∑

ω2el

Задавая вектор k, из этого уравнения можно опре-

ε_L (ω) = 1 +

(18)

ω20el - (ω² - iωω_cl)

делить волновое число:

l=1

Здесь квадрат частоты ω2el

= e²N_l/ (ε₀m_l) опре-

k²⁰ =

деляется концентрацией N_l и массой заряженных

(17)

частиц сорта l, а частота ω_cl релаксации импульса

∑

4δ²b (ε - 1)

|G_klm (k)|²

(частота столкновений) определяет уширение спект-

a³

+k2yl +

-k²

ральной линии. Квантовое описание на основе по-

k,l,m=-∞

луклассического приближения дает такую же фор-

Уравнение (17) можно решать итерационно. Оно

му ДП [6, 13, 32]:

также подходит для металлической рамки, если в

качестве ε взять ДП металла. В случае наличия ос-

e²

ε_L (ω) = 1 +

новы ДП ε необходимо сделать замены k²⁰ → k²⁰ ε,

ε₀m_e

∑∑

∑

ε → ε/ε. Считая k = x₀k_x, k_x ≈ k₀ ≪ π/a, учи-

f^νlk

N^ν

(19)

тывая только большой нулевой член суммы, полу-

⌊

⌋

ν k

l=k

(ων0lk)² - (ω² - iωω^νlk)

чим k²⁰

≈ k2x/

1 + 4δ²b(ε - 1)/a³

. Приближение

справедливо, если 4δ²b (ε - 1) /a³ ≪ 1. В этом слу-

Здесь суммирование по индексам ν идет по всем сор-

чае k_x ≈ k₀. Конечно, уравнение (17) весьма при-

там атомов, по индексам k — по всем энергетиче-

ближенное. Задавая поле в виде циклического то-

ским уровням, по l — по всем разрешенным пере-

ка поляризации, мы наложили условия симметрии:

ходам с уровня k. Соответственно N^νk — число ато-

E_x (x, y, z) — четная функция по x, z и нечетная

мов типа ν в состоянии k, ω0ν = ω^νlk = (E_k - E_l) ℏ,

по y, а E_y (x, y, z) — четная функция по y, z и нечет-

ω^νlk — ширина спектральной линии перехода k → l,

ная по x. В силу симметрии такое решение может

связанная со временем жизни. Если в возбужден-

существовать, но могут быть решения с другими со-

ных состояниях атомов больше, чем в основном,

отношениями четности-нечетности. Решая задачу в

то ε^′′ < 0 и среда является активной. Сила осцилля-

общем виде, например, методом конечных элемен-

тора f^νlk = 2m_e |ων0lk| |d^νlk|² / (3e₂ℏg_l) дипольного пе-

тов, мы можем найти волны всех возможных сим-

рехода для атома с дипольным моментом d^ν опре-

метрий и соответствующие им дисперсии. Уравне-

деляется квадратом модуля матричного элемента

ние (17) при полностью четном поле будет четным.

dνlk = -e 〈ψl |r| ψk〉 [6, 32]. В общем случае при вы-

При дифракции волн на метаматериале конечных

числении (19) следует суммировать по всем элект-

209

ЖЭТФ, вып. 2

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

рическим мультипольным переходам, удовлетворя-

ющим правилам отбора. Нетрудно видеть, что мо-

жет выполняться условие ε^′L (ω) < 0. Однако (19)

не дает правильное значение ДП, поскольку полу-

чено при воздействии поля плоской волны в ва-

кууме на отдельный атом и путем суммирования

вкладов атомов в поляризацию. Необходим учет

среднего поля [6]. Часто в работах формулу ти-

па (18) пишут для μ_e, пытаясь получить μ^′e < 0

из решений задач дифракции, используя пакеты

компьютерного моделирования. Но формула типа

(18) не имеет отношения к μ_e. В квантовом рас-

смотрении, как известно, можно представить воз-

мущенный гамильтониан в виде

Ĥ = Ĥ₀ + Ĥ^′, где

Ĥ′

= -(e/m) p · A^e + (eA^e)2 /(2m) +

Ĥ_s

— свя-

занное с полем B = ∇ × A возмущение [32, 33].

Здесь

Ĥ_s

= μ_Bg_sB · ŝ/μ₀, A^e = μ₀A, A — ранее

введенный вектор-потенциал, gs ≈ 2 и для прос-

Рис. 5. Дисперсия в метаматериале из прямоугольных ди-

тоты рассмотрен один электрон в атоме со спи-

электрических рамок

ном ŝ. Пренебрегая вторым членом и учитывая

свойства коммутации, можно записать [33]

Ĥ′

(

)

=μ_B

^L+2Ŝ /μ0, где μB — магнетон Бора,

^Ŝи

^L—

ных частот. При комнатной температуре длина

полные операторы спинового и орбитального момен-

тов атома (включая ядро). Магнитная восприимчи-

свободного пробега электронов в металле порядка

десятков нанометров. Кроме того, в отличие от

вость зависит от полного магнитного момента ато-

ма, и если он ненулевой, слагается из поляризацион-

спина в атоме магнитный момент контура не может

менять направление ориентации в зависимости от

ного парамагнетизма, прецессионного диамагнетиз-

ма и ориентационного парамагнетизма [34]. В сла-

изменения угла магнитного поля. Такой наведенный

бых полях для определения вклада атома в маг-

магнетизм на низких частотах является диамагне-

тизмом [26]. Можно получить формулу типа (18)

нитную поляризацию слеует выислять матричные



без резонансов для μ-1e, поскольку на движущиеся

элементы m_lk = -μ_B ψ_l

2^Ŝ

+L ψk . Ротационные

заряды действует поле B (ω). Вклад в поляри-

спектры полярных молекул в газе и жидкости также

зацию от МА пропорционален B (ω), а не H (ω).

могут формировать магнитный момент. Однако это

Классический подход легко можно применить к

обычно микроволновый диапазон с сильным уши-

интегральному уравнению из проволочных мик-

рением линий, и какого-либо вклада в МП не воз-

роколечек (рис. 4в). Удобно расположить их на

никает. В оптическом диапазоне природные веще-

ребрах кубического кристалла. Часто пытают-

ства не демонстрируют существенное отличие МП

ся рассматривать резонансные структуры типа

от единицы [9, 34]. ИС из металлических наноча-

разомкнутых микроколечек с емкостным зазором.

стиц, представляющих собой квантовые точки или

У такого RCL-контура есть резонансная часто-

трехмерные квантовые ящики, целесообразно моде-

та (LC)-1/2 вынужденного резонанса, но она очень

лировать квантовыми методами. Определяя для них

высокая, поскольку краевая емкость мала. Для

уровни энергии и волновые функции или исполь-

снижения ω₀ следует повышать емкость. Кольцо

зуя теорию функционала плотности, можно вычис-

с емкостью является магнитным и электрическим

лить плотность тока, возникающего под действием

диполями, ориентированными нормально и каса-

поля B (ω). Это поле также действует и на классиче-

тельно плоскости кольца. Электрический диполь

скую рамку с током, поэтому в слабых полях P^m (ω)

нормален зазору, что делает кольцо несиммет-

пропорционально B (ω), а не H (ω).

ричным излучателем в его плоскости. Поэтому

Рассматривая классические контуры с токами,

кристалл (рис. 4в) с разомкнутыми кольцами не

возбуждаемыми магнитной индукцией B (ω)

изотропный. Не спасает и двойное разомкнутое

= μ₀μ_eH (ω), следует учесть, что ток в контуре

кольцо. Необходимы четверные разомкнутые коль-

в отсутствие поля затухает и не имеет резонанс-

ца с поворотом зазоров на

45 %, размещенные

210

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

на ребрах кубического фотонного кристалла. В

7. УРАВНЕНИЕ ФРЕНЕЛЯ И ЧИСЛЕННЫЕ

силу правила Ленца такая ИС на низких частотах

РЕЗУЛЬТАТЫ

обладает диамагнетизмом. Важно синхронизиро-

Периодические фотонные кристаллы описыва-

вать углы поворота всех колец. Реально такой

ются уравнением Френеля. Это дисперсионное урав-

3D-фотонный кристалл возможен в радиодиапа-

нение, построенное по результатам гомогенизации.

зоне, где можно использовать и сосредоточенные

В классическом случае имеют место два подхода

емкости. Создание такой ИС в микроволновом

к описанию эффективных параметров ИС: симмет-

диапазоне уже проблематично. Вблизи резонансной

⌊

(

)⌋

ричный в рамках полей E, H и индукций D, B, а

частоты μ-1e (ω)

1 + ω2m/

ω²⁰ -

ω² - iωω_c

также несимметричный с использованием трех век-

Рассмотрим μ^′e. Имеем

торов E, D, B [8,9]. В этом разделе будем придер-

(

)

ω2m

ω²

-ω² +ω²

живаться первого подхода. В нем в общем случае

μ^′e = 1 -

бианизотропной ИС вводятся четыре тензорных па-

(ω²⁰ - ω² + ω2m)² + (ωω_c)²

раметра: тензоры ДП ε(ω, k) и МП μ (ω, k) и тен-

(

)₂

зоры кросс-поляризаций

ξ (ω, k),

^ζ (ω, k). При рас-

Условие μ^′e < 0 имеет вид

ω2m - ω2c

> 4(ω_mω_c)².

(

)₂

пространении плоских волн E = E₀ exp (iωt - ik · r)

Обозначим Δ =

ω2m - ω2c

> 4(ω₀ω_c)². При слабой

(или при комбинации поля из таких волн) и при сим-

диссипации ω2m - ω2c > 2ω₀ω_c или ω_m >

^√2ω₀ω_c. В

метричном рассмотрении имеют место связи

этом случае имеется область отрицательных μ^′e < 0,

расположенная при

D(ω, k) = ε₀ ε(ω, k)E (ω, k) + c-1

ξ (ω, k) H (ω, k) ,

ω2m - ω2c

B (ω, k) = μ₀ μ (ω, k) H (ω, k) + c-1

ζ (ω, k) E (ω, k) ,

ω²⁰ +

-Δ<ω² <ω²⁰ +

+ Δ.

а также уравнения Френеля вида [1, 9]

Обозначая Ω = ω² - ω²⁰, имеем уравнение Ω² -

⌊(

)

(

)

⌋

(

)

-Ω

ω2m - ω2c

+ (ω₀ω_c)² = 0 для определения гра-

det k-10k+

μ-1

k-10k - ζ + ε = 0,

ничных частот, на которых μ^′e = 0. С учетом дисси-

⌊(

)

(

)

⌋

(20)

det k-10k-ζ

ε-¹

k-10k +

+μ .

пации таких частот

⌊

⌋

(

)₂

частоте μ^′e = 1 - 1/ 1+

ω₀ω_c/ω2m

> 0. В этом

Здесь ∇×E = -ik×E =

^kE, а матрица

k определяет

случае при ω₀ω_c = ω2m получаем μ^′e = 1/2. С рос-

оператор ротор:

том диссипации μ^′e при резонансе растет и стремит-

⎡

⎤

ся к 1, а при стремлении диссипации к нулю также

-k_z k_y

стремится к нулю. При низких частотах имеем 0 <

⎥

k₌ -i^⎣ k_z

-k_x⎦.

< μ^′e < 1, т.е. диамагнетизм, и μ^′e → 1 при ω → 0.

−k_y k_x

На очень больших частотах (ω² ≫ ω²⁰ + ω2m) име-

(

)

ем μ^′e = 1 +

ω2m + 2ω2c

/ω² → 1, т. е. реализуется

Если рассматривать изотропные метаматериалы в

парамагнетизм. Отметим, что мы ввели резонанс-

смысле работы [4], т. е. считать все эффективные

ную частоту искусственно на основе теории цепей

материальные параметры скалярами, то ИС стано-

и воспользовались формулой, которую еще надо по-

вится в общем случае биизотропной со скалярны-

лучить, используя электродинамику. То есть надо

ми па⌊аметрами ε, μ, ξ, ζ ⌋и уравнением Френе-

строго решать задачу для разомкнутых колец.

ля det

k² +^k (ξ - ζ) +

Ik²⁰εμ

= 0. Если же кросс-

Таким образом, отрицательная МП может быть

реализована в резонансных ИС, но только в резо-

поляризацией можно пренебречь, а анизотропией

нансной области, где существенна ПД, и описание

нельзя, то уравнение Френеля приобретает вид

двумя скалярными величинами ε_e и μ_e не является

⌊

⌋

полностью корректным. В низкочастотной области,

det

ε-¹kμ-1k -

Ik20

= 0,

⌊

⌋

(21)

где можно ввести эти параметры, ИС ведет себя как

det

μ-1kε-1 -

Ik²

= 0.

диамагнетик: 0 < μ_e < 1. Поскольку размер части-

цы существенно меньше a, указанные резонансы ле-

В полностью изотропном случае (ε = ε_ef

I,μ=μ_efI)

жат обычно выше первого брэгговского резонанса,

⌊

⌋

где уже изотропное приближение неприменимо. На-

уравнение Френеля det

k² -

Ik²⁰εμ

= 0 экви-

веденный магнетизм, однако, вполне имеет место в

валентно простейшему дисперсионному уравнению

широком диапазоне частот.

для обыкновенной волны k² = k²⁰εμ.

211

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

Для ИС из шариков двух типов рис. 4г оценим

резонансные частоты, что, конечно, не означает до-

стижимость однородного метаматериала, описыва-

емого ДП и МП. Для обоснования этого следует

задать в шариках меньшего размера только поле

E_ϕ = const и расположить эти шарики на ребрах.

Уже отмечалось, что такого возбуждения при ди-

фракции волны на конечном образце может не быть.

Поскольку ориентация осей различная, шарики да-

ют вклад в магнитный момент по трем направле-

ниям. Такие шарики подобны рассмотренным коль-

цам и создают магнитный момент. Конечно, это

приближение: поле в шариках следует разлагать по

всем сферическим гармоникам [18,25]. Пусть шари-

ки большего размера расположены в углах. Эти ша-

рики создают электрический момент. ДП шариков

считаем одинаковой. Для разложения поля можно

Рис. 6. Эффективная МП метаматериала из прямоуголь-

взять несколько сферических гармоник относитель-

ных диэлектрических рамок по формуле (15), соответству-

но трех осей. В такую ИС мы искусственно при-

ющая дисперсии рис. 5

внесли симметрию. Далее следует рассчитать воз-

буждаемое шариками поле с учетом периодической

(

)

ФГ (1) и решить интегральное уравнение. Такое ре-

ля Бруггемана. Поскольку C = (4π/3)

r³¹ + r³²

/a³,

шение учитывает взаимовлияние всех шариков. За-

при r₂ = 2r₁ имеем r₁ ∼ a·10-4. При ε_e = 11 должно

тем следует вычислить электрический и магнитный

быть Λ < 1 мм, т. е. по крайней мере a < 0.1 мм, что

моменты и усредненные поля. В природных средах

дает оценку r₁ ∼ 10-5 мм. Минимальное расстоя-

√

размер атомов r ≈ 0.05 нм, и при длине световой

ние между рассеивателями a

3/2, поэтому взять r₁

волны λ ≈ 500 нм и ДП основы ε ≈ 10 и менее при-

больше 0.01 мм без существенного нарушения раз-

ходим к отношению размера к длине волны в сре-

реженности нельзя. Однако такое увеличение нару-

де r/Λ ≈ 3 · 10-4. Для кристаллов это отношение

шает условие ε_e ≈ ε и уменьшает Λ. Нетрудно ви-

имеет такой же порядок. На атомы действует ло-

деть на основе точных формул [36], что резонанс-

кальное поле, а прецессия спинов дает магнитный

ные частоты таких изолированных наношариков не

момент, в линейном приближении пропорциональ-

лежат в заявленном авторами идеи диапазоне. В

ный компоненте локального магнитного поля. Для

свободном пространстве основной (низкочастотной)

возбуждения рамки с током или шарика важен их

является магнитная мода H01δ, при этом индексы

размер, сравнимый с длиной волны. Магнитный мо-

(кроме азимутального) нельзя считать целыми (в

мент не пропорционален локальному полю, а намаг-

силу радиационных потерь их даже нужно считать

ниченность зависит от k. Для фотонного кристалла

комплексными). В фотонном кристалле нет радиа-

с a ≈ 50 нм уже a/Λ ≈ 0.1, при этом для разре-

ционных потерь, но есть диполь-дипольное взаимо-

женности необходимо, по крайней мере, использо-

действие рассеивателей, которое необходимо строго

вать частицы с r < 5 нм. В области 100 ГГц при

учитывать.

длине волны 3 мм, ε = 10 и ε = 400 оценим ми-

Низкочастотная дисперсионная ветвь, соответ-

нимальную Λ. Поскольку изотропность предполага-

ствующая уравнению (17), приведена на рис. 5. На

ет работу вдали от брэгговских резонансов, хоро-

рис. 6 (кривая 1) даны результаты вычисления МП

шим и точным подходом для этого являются мето-

по формуле (15). В суммах формул (15) и (17) име-

ды для хаотических метаматериалов типа смесей:

ют место полюсы, которые определяются условием

формулы смешения для моделей эффективной сре-

+k2yl +k2zm = k²⁰. Они приводят к нулям в правых

ды, метод компактных групп, методы теории перко-

частях (15) и (17). Кроме того, имеют место нули

ляции (протекания) и т. п. [7, 35]. Из формулы Гар-

у функций F_mlk (k) и |G_mlk (k)|². Они реализуются

нетта (ε_e - ε) / (ε_e + 2ε) = C (ε - ε) / (ε + 2ε) имеем

при достаточно больших |k| и приводят к полюсам.

оценку для концентрации C диэлектрических шари-

Такие резонансы в МП нефизические и происходят

ков при условии ε_e = ε + 1 = 11: C = 1/30. Близ-

при достижении запрещенной зоны, т. е. там, где мо-

кий результат дает и формула эффективного по-

дель неприменима. Избежать полюсов можно, учтя

212

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

потери или построив более строгую модель. В низ-

кочастотной области |k| ≈ k₀ нуль компенсирует по-

люс и МП практически равна единице, а модель де-

монстрирует бесконечно малый диамагнетизм. От-

рицательные значения не достигаются. Расчеты в

этой области выполнены при введении бесконеч-

но малого тангенса диэлектрических потерь 10-10.

При больших |k| для вычисления МП использова-

лись связи kx = kx (k0) и k = kx + ky = k (k0),

определяемые из уравнения (17), и брался тангенс

потерь 10-2. Заметим, что можно построить дис-

персию МП по результатам гомогенизации на ос-

нове уравнения Френеля μ = (k_x/k₀)². Такая мо-

дель при |k| ≈ k₀ дает слабый парамагнетизм. В

точке X имеем μ = 1.061. При переходе на обрат-

ную ветвь k_x → π/a - k_x замедление и МП силь-

но растут. В запрещенной зоне компонента k_x яв-

ляется мнимой, а МП — отрицательной. При пере-

ходе через запрещенную зону на высокочастотную

(оптическую) дисперсионную ветвь при изменении

Рис. 7. Дисперсия (низшие частотные ветви) в кубическом

только k_x эта МП становится сначала отрицатель-

фотонном кристалле с кубическими МА a0 × a0 × a0 в ди-

ной, затем положительной, но имеет значение мень-

электрической основе ε = 3.0. ДП включений ε = 6.0 для

ше единицы. Как для рис. 6 в резонансных областях,

кривых 1-4, ε = 0.5 (5), ε = -1.0 (6), ε = 3.1 (7). Кривым

так и здесь отличные от единицы резонансные зна-

1, 4-7 соответствует a0 = 0.1a, a0 = 0.2a (2), a0 = 0.5a (3)

чения МП не имеют физического смысла, хотя эти

значения и описывают дисперсию моделью Френе-

ля. При дифракции на образце фотонного кристал-

ла для использования такой модели важно, в какой

точке дисперсионной поверхности при данной часто-

те находится вектор k с учетом всех волн, включая

и отраженные, что само по себе сложно.

На низких частотах можно построить другие за-

висимости замедления на основе эффективной ди-

электрической среды в виде смеси. Для нее имеем

коэффициент заполнения C = 16bδ²/a³ и коэффи-

циенты деполяризации (см. [7, 37]) L_x = L_y = 1/2 -

- δ/b, L_z = 2δ/b. Используя формулу Гарнетта, для

тонкой диэлектрической рамки получим [1]

C (ε - 1)

ε_exx = ε_eyy = 1 +

≈

1 + Lx (ε - 1 - C)

16bδ² (ε - 1)/a

≈1+

1 + (ε - 1 - 16bδ²/a³)/2

Рис. 8. Гомогенизированные значения эффективной ДП

16bδ² (ε - 1)/a

ε_ef для конфигураций (1-7), соответствующих рис. 7. Кри-

ε_ezz = 1 +

≈

1 + 2δ(ε - 1 - 16bδ²/a³)/b

вая 8 соответствует ε = -1.0, a0 = 0.3a, ε = 1.0

≈ 1 + 16bδ2 (ε - 1)/a³.

Результаты вычисления низкочастотных диспер-

Таким образом, на низких частотах ИС описывает-

сионных ветвей для (17) с эффективной ДП для

ся диэлектрическими свойствами, а волна с k_z = 0

кубического кристалла с кубическими диэлектри-

идет почти со скоростью света. Теперь уравнение

ческими МА даны на рис. 7 и 8. Использовалась

Френеля при k_y = k_z = 0 имеет вид k_x = k₀

√ε_ezz.

приближенная модель в виде аппроксимации ком-

213

М. В. Давидович

ЖЭТФ, том 159, вып. 2, 2021

понент E_x, E_y, E_z кусочно-постоянными элемента-

что большая МП не имеет обычного физического

ми с симметричным распределением поля, поэтому

смысла. В низкочастотном пределе при малом фа-

уравнение решалось в первом октанте куба. Исполь-

зовом сдвиге на ячейку МП принимает локальный

зовалось по две и три аппроксимации на компонен-

смысл.

ту, поэтому размерность задачи равнялась либо 8,

Для поддерживающего волну тока поляризации

либо 27. Гомогенизация выполнялось по формуле,

J_p (даже замкнутого) намагниченность можно не

аналогичной (15). Заметим, что простой алгоритм

вводить, а ввести только электрическую поляриза-

можно получить, взяв аппроксимации полей типа

цию 〈P^e〉 = ε₀ (ε - 1) 〈E〉_V . Усреднение здесь выпол-

E_x (x, y, z) = E0x cos(α_xx)cos(β_xy)cos(γ_xz). Тогда

нено по МА. Такой дипольный момент отличается

дисперсионное уравнение получается в виде функ-

от (5) и является более приближенным. В низкоча-

ционала, зависящего кроме k₀ и k еще от 12 парамет-

стотном пределе поле постоянно и 〈E〉_V = V 〈E〉/a³.

ров. Их можно сократить до 8, наложив условие со-

Формально можно вводить разные модели ИС, ис-

леноидальности и условия удовлетворения каждой

пользуя тензорные или скалярные материальные

из компонент волновому уравнению. Однако такая

параметры, а затем определять их из условий гомо-

задача является нелинейной.

генизации. В качестве таковых могут фигурировать

свойства волн в метаматериале, параметры рассея-

ния, поляризуемости. Неоднозначность обратной за-

8. ВЫВОДЫ

дачи гомогенизации требует решения переопреде-

Рассмотрены структуры периодических метама-

ленной системы уравнений. В частных случаях чис-

териалов, которым могут соответствовать изотроп-

ло уравнений и неизвестных может совпадать. Опре-

ные эффективные материальные параметры, в том

деление МП такими способами не гарантирует, что

числе и с двумя проницаемостями ε_e и μ_e. Посколь-

имеет место локальная связь магнитной индукции

ку все МА (исключая ферромагнитные частицы во

и магнитного поля. Формально введение магнит-

внешнем магнитном поле) можно описать ДП, сим-

ных свойств основано на том, что в низкочастотном

метричный кубический кристалл также моделиру-

пределе рассеяние плоской волны на малой частице

ется только диагональным тензором с одинаковыми

можно описать как возникновение в ней электриче-

компонентами ε_e (k₀, k) [1,6]. Поэтому введение μ_e —

ского и магнитного диполей [7, 17, 37]. Однако для

искусственный прием. Исключение могут составить

выполнения граничных условий магнитный диполь

магнитные МА из ферритов и магнитных металлов,

можно не вводить, а ввести еще один ортогональный

магнетизм которых обусловлен прецессией неком-

электрический диполь. Мы не рассмотрели ряд важ-

пенсированных магнитных моментов атомов или до-

ных вопросов: влияние квадрупольных и высших

менов во внешнем магнитном поле. Такие ИС могут

мультипольных моментов на поляризацию, плотно

демонстрировать резонансные свойства и изменение

упакованные (не разреженные) фотонные кристал-

знака компонент МП, но они анизотропны и обычно

лы, сильно нелокальные ИС и ряд других.

проявляют такие свойства в микроволновом диапа-

Итак, однородные периодические ИС, описывае-

зоне и ниже. Там имеет место гиротропия.

мые ДП, возможны в низкочастотном пределе при

Так называемый наведенный магнетизм демон-

пренебрежении ПД. При описании их ДП и МП по-

стрирует зависимость намагниченности от k даже

следняя близка к единице, особенно в оптическом

в низкочастотном пределе. При этом μ_e (0, 0) = 1.

диапазоне. Однородные метаматериалы с отрица-

Показано, что для кольцевых ИС в низкочастот-

тельной ДП также возможны. Это металлы при

ном пределе |k| ≈ k₀ (0, 0) ≪ π/a имеет место сла-

сверхнизких температурах на частотах вплоть до

бый диамагнетизм (см. [26]). Рассмотрим это качест-

плазменных. Смеси металлических и диэлектриче-

венно на примере одиночной проволочной рамки с

ских МА в низкочастотном пределе также могут

площадью S. Имеем A =

^Ĝ(J), H = ∇ × A =

демонстрировать отрицательную ДП. Наведенный

= -jk ×

^Ĝ (J), dM = r × Jdl. Для плоской круг-

магнетизм возможен, но соответствующие ему мета-

лой рамки c током I = |J| δ² получаем M = n₀SI.

материалы не являются изотропными и могут быть

Диамагнетизм обеспечивает правило Ленца. Усред-

описаны другими материальными параметрами. В

ненное магнитное поле не является локальным, по-

оптике МП метаматериалов практически равна еди-

этому нельзя ввести локальные материальные па-

нице, что согласуется с аналогичным выводом для

раметры. Большие значения МП связаны с тем, что

природных веществ [7, 30]. Для периодических ИС

〈H_z (k₀, k)〉 ≈ 0 для некоторых точек дисперсионной

на поставленный в заглавии вопрос следует ответить

поверхности. Нелокальное магнитное поле означает,

отрицательно.

214

ЖЭТФ, том 159, вып. 2, 2021

Возможны ли изотропные метаматериалы. . .

Финансирование. Работа осуществлена при

19.

С. М. Рытов, ЖЭТФ 29, 605 (1955).

поддержке Министерства образования и науки Рос-

20.

T. Gric, Waves in Random and Complex Media 33(4)

сии в рамках выполнения государственного задания

(2019),

https://doi.org/10.1080/17455030.2019.

(проект № FSRR-2020-0004).

1656846.

21.

М. В. Давидович, Известия Сарат. унив. Серия

ЛИТЕРАТУРА

Физика 11(1), 42 (2011).

М. В. Давидович, УНФ 189, 1250 (2019).

22.

М. В. Давидович, КЭ 47, 567 (2017).

А. П. Виноградов, А. В. Дорофеенко, С. Зухди,

23.

А. В. Вашковский, Э. Г. Локк, УФН 176, 557

УНФ 178, 511 (2008).

(2006).

М. Я. Сушко, С. К. Криськив, ЖТФ 79(3), 97

24.

В. П. Макаров, А. А. Рухадзе, УФН 189, 519

(2009).

(2019).

К. Р. Симовский, Опт. и спектр. 107, 766 (2009).

25.

Г. Т. Марков, А. Ф. Чаплин, Возбуждение элект-

ромагнитных волн, Радио и связь, Москва (1983).

М. А. Ремнев, В. В. Климов, УНФ 188, 169 (2018).

26.

М. В. Давидович, Письма в ЖЭТФ 108, 299

В. М. Агранович, Ю. Н Гартштейн, УНФ 176, 1051

(2018).

(2006).

27.

И. Б. Вендик, О. Г. Вендик, М. С. Гашинова, Пись-

А. П. Виноградов, Электродинамика компо-

ма в ЖТФ 32(10), 30 (2006).

зитных материалов, Эдиториал УРСС, Москва

(2001).

28.

O. Vendik, I. Vendik, I. Kolmakov, and M. Odit,

Opto-Electron. Rev. 14(3), 179 (2006).

А. П. Виноградов, УФН 172, 363 (2002).

29.

L. Jylha, I. Kolmakov, S. Maslovski, and S. Tretya-

Л. Д. Ландау, Е. М. Лифшиц, Теоретическая фи-

kov, J. Appl. Phys. 99, 043102 (2006).

зика. Электродинамика сплошных сред, Наука,

Москва (1982).

30.

A. Ahmadi and H. Mosallaei, Phys. Rev. A 77,

045104 (2008).

10.

R. Merlin, Proc. Nat. Acad. Sci. 106, 1693 (2009).

31.

A. Krasnok, S. Makarov, M. Petrov, R. Savelev,

11.

A. H. Aly, M. Ismaeel, and E. Abdel-Rahman, Opt.

P. Belov, and Yu. Kivshar, Proc. SPIE 9502, 950203

Photon. J. 2, 105 (2012).

(2015).

12.

М. В. Давидович, И. А. Корнев, Компьютерная оп-

32.

М. О. Скалли, М. С. Зубайри, Квантовая опти-

тика 43, 765 (2019).

ка, под ред. В. В. Самарцева, Физматлит, Москва

13.

А. И. Ахиезер, И. А. Ахиезер, Электромагне-

(2003).

тизм и электромагнитные волны, Высшая шко-

33.

Л. Д. Ландау, Е. М. Лифшиц, Квантовая ме-

ла, Москва (1985).

ханика (нерелятивистская теория), Физматлит,

14.

R. D. Graglia, P. L. E. Uslenghi, and R. E. Zich, IEEE

Москва (1963).

Trans. AP-39(1), 83 (1991).

34.

Я. Г. Дорфман, Магнитные свойства и строение

15.

E. O. Kamenetskii, Phys. Rev. E 57, 3563 (1998).

вещества, Гостехиздат, Москва (1955).

16.

E. O. Kamenetskii, Microw. Opt. Technol. Lett. 19,

35.

А. А. Снарский, УФН 177, 1341 (2007).

412 (1998).

36.

P. Guillon and Y. Garault, IEEE Trans.

17.

L. Lewin, Proc. Inst. Elec. Eng. 94(27), 65 (1947).

MTT-25(11), 916 (1977).

18.

Л. Левин, Современная теория волноводов, Изд-во

37.

R. Е. Соllin, Field Theory of Guided Waves, IEEE

иностр. лит., Москва (1954).

Press, New York (1991).

215