ЖЭТФ, 2022, том 162, вып. 1 (7), стр. 14-26

ОПРЕДЕЛЕНИЕ ДИСПЕРСИОННЫХ ХАРАКТЕРИСТИК

ИНТЕГРАЛЬНЫХ ОПТИЧЕСКИХ МИКРОРЕЗОНАТОРОВ

И ГЕНЕРАЦИЯ В НИХ КОГЕРЕНТНЫХ ОПТИЧЕСКИХ

ЧАСТОТНЫХ ГРЕБЕНОК

Н. Ю. Дмитриев^a,b*, А. С. Волошин^c, Н. М. Кондратьев^a, В. Е. Лобанов^a,

a,d

К. Н. Миньков^a, А. Е. Шитиков^a, А. Н. Данилин^a,d, Е. А. Лоншаков^a, И. А. Биленко

^a Российский квантовый центр

121205, Москва, Россия

^b Московский физико-технический институт

141701, Долгопрудный, Московская обл., Россия

^c Swiss Federal Institute of Technology Lausanne (EPFL)

1015, Lausanne, Switzerland

^d Физический факультет, Московский государственный университет им. М. В. Ломоносова

119991, Москва, Россия

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

после переработки 17 февраля 2022 г.

Принята к публикации 17 февраля 2022 г.

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

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

рентных частотных гребенок служат высокодобротные микрорезонаторы. Ключевыми характеристиками,

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

резонатора. В данной работе разработана и верифицирована оригинальная методика, позволяющая с

высокой точностью определить величину дисперсионных коэффициентов как в кристаллических, так и в

интегральных микрорезонаторах с межмодовым интервалом от гигагерц до терагерц. Продемонстриро-

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

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

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

DOI: 10.31857/S0044451022070021

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

EDN: DRLYOD

ческая частотная гребенка, представляющая собой

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

линий. Оптические частотные гребенки в высоко-

1. ВВЕДЕНИЕ

добротных оптических микрорезонаторах впервые

были продемонстрированы в 2007 году [12]. Генера-

Высокодобротные оптические микрорезонаторы

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

с модами шепчущей галереи, впервые продемон-

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

стрированные более 30 лет назад [1], открыли уни-

ми, происходит за счет процесса нелинейного четы-

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

рехволнового взаимодействия [13-16]. Процесс гене-

ния широкого круга оптических явлений и до сих

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

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

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

[2-11]. Одним из наиболее активно развивающих-

объемных кристаллических и интегральных [17-19].

ся направлений является генерация в них такого

Хотя оптические частотные гребенки известны и

^* E-mail: nkdmitriev@gmail.com

успешно применяются в науке и технике уже более

ЖЭТФ, том 162, вып. 1 (7), 2022

Определение дисперсионных характеристик. . .

hν1

hν3

hν2

hν4

(

) *'

!"#$

)+,+

!"#$

Рис. 1. Схематичное представление процесса генерации частотной гребенки: в фотонный чип со сверхдобротным мик-

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

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

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

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

исходит каскадный процесс генерации боковых линий, отстоящих друг от друга строго на величину области свободной

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

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

гребенка), во временном представлении имеет вид последовательности сверхкоротких импульсов (солитонов)

20 лет [20,21], активное развитие интегральной фо-

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

тоники открыло возможности существенной мини-

ющегося в среде сигнала, и дисперсией групповых

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

скоростей (ДГС), с одной стороны, и между потеря-

волило существенно расширить область их приме-

ми и накачкой, с другой. Нелинейность третьего по-

нения. На сегодняшний день микрогребенки успеш-

рядка, обычно называемая керровской, в основном

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

зависит от природы используемого материала и его

зи [22, 23], оптические сверхбыстрые дальномеры

свойств. ДГС же зависит одновременно и от свойств

[24,25], спектроскопия в режиме реального времени

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

[26], калибровочные эталоны [27, 28], атомные часы

резонатора.

[29], интегральные источники микроволнового излу-

В зависимости от знака коэффициента диспер-

чения [30, 31] и т. п.

сии второго порядка выделяют два типа ДГС, нор-

Особый интерес представляют когерентные час-

мальную и аномальную. Динамика нелинейных про-

тотные гребенки, соответствующие сверхкоротким

цессов в микрорезонаторе зависит от типа ДГС

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

на длине волны накачки [34]. Для генерации со-

нии. Такие сигналы характеризуются наименьшим

литонных частотных гребенок, как правило, тре-

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

буется аномальная ДГС. Обладая информацией о

ральной огибающей. Генерацию таких частотных

ее величине, можно однозначно определить пара-

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

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

можно наблюдать при накачке высокодобротного

ных значений мощности накачки и отстройки час-

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

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

ем при выполнении определенных условий [32, 33]

тора [32]. Более того, такой параметр, как ДГС, яв-

(рис. 1). Данное явление становится возможно при

ляется одним из важнейших параметров в физике

соблюдении баланса между нелинейностью третьего

сверхкоротких импульсов и играет ключевую роль в

порядка, которая выражается зависимостью показа-

большинстве нелинейных процессов. В связи с этим

Н. Ю. Дмитриев, А. С. Волошин, Н. М. Кондратьев и др.

ЖЭТФ, том 162, вып. 1 (7), 2022

! "

,'6$

245'.$&'2$

789!:

f-2

f-1

#$%&'&$()*+)

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

ду ОСД микрорезонатора (D1/2π, штриховые линии) и межмодовым интервалом при учете вклада дисперсии второго

порядка (см. пунктирные линии) в случае аномальной дисперсии D2/2π > 0

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

мым микрорезонатором [2, 11].

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

характеристик систем, а также иметь возможность

2. МЕТОДИКА ЭКСПЕРИМЕНТА И

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

ИЗМЕРИТЕЛЬНЫЙ СТЕНД

резонатора [35, 36]. В ряде случаев прямое измере-

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

2.1. Экспериментальная оценка

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

дисперсионных характеристик

рорезонаторов с большой величиной области свобод-

Общая идея экспериментальной оценки диспер-

ной дисперсии, превышающей сотню гигагерц.

сионных характеристик микрорезонатора заключа-

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

ется в определении положения всех его собственных

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

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

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

щихся к одному пространственному семейству, и в

микрорезонаторов и проведена его верификация.

последующем вычислении дисперсионных коэффи-

Преимущества рассматриваемого метода заключа-

циентов D_i. Вычисление дисперсионных коэффици-

ются в его гибкости и доступности. Рассмотренный

ентов осуществляется путем аппроксимации полу-

метод позволяет проводить измерение ДГС прак-

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

тически для любого диапазона длин волн. Поми-

формулой зависимости частоты моды от ее поряд-

мо этого, он подходит для измерения дисперсион-

кового номера µ [32, 35]:

ных характеристик как интегральных микрорезо-

D₂

D₃

наторов с большим значением области свободной

ω_µ = ω₀ + D₁µ +

µ² +

µ³ +

дисперсии (ОСД), так и кристаллических с рекорд-

D₄

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

µ⁴ + . . . ≡ ω₀ + D₁µ + D_int,

(1)

зонансами. Вдобавок рассмотренная методика ка-

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

где ω₀

собственная частота моды накачки,

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

D₁/2π

ОСД резонатора, D₂/2π

показатель

пример, для оценки добротности. Также в работе

дисперсии второго порядка, определяющий ДГС,

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

D₃/2π, D₄/2π, . . . коэффициенты, относящиеся к

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

более высоким порядкам дисперсии. Будем предпо-

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

лагать, что D₁ ≫ D₂ ≫ D₃ ≫ . . ., и ограничимся

характеристики полностью соответствуют рассчи-

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

танным на основе измеренных с помощью разрабо-

но часто под ОСД подразумевают интервал между

танной методики дисперсионных параметров микро-

соседними модами резонатора, азимутальные числа

резонатора. Для генерации частотной гребенки в ка-

которых отличаются на 1, однако, как видно из со-

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

отношения (1), интервал между соседними модами

ный лазерный диод, стабилизированный тем же са-

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

ЖЭТФ, том 162, вып. 1 (7), 2022

Определение дисперсионных характеристик. . .

В данной работе межмодовым будем называть

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

интервал, который складывается из ОСД D₁/2π и

остро стоит вопрос высокоточной калибровки в про-

добавки D₂/2π, обусловленной наличием дисперсии

цессе проведения измерений. Данная задача может

второго порядка (см. рис. 2). В случае резонатора

быть решена путем использования прецизионного

с нулевой дисперсией (D_i = 0, i ≥ 2) межмодовый

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

интервал и ОСД совпадают: все моды будут рас-

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

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

лонов, таких как источники гребенок [40]. Однако

равном ОСД. Аномальной ДГС соответствует

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

случай D₂ > 0, а нормальной случай D₂ < 0.

статков. Во-первых, подобное оборудование, облада-

Для генерации солитонных частотных гребенок,

ющее необходимым уровнем точности, стоит крайне

как правило, требуется аномальная ДГС. В ряде

дорого, и многие лаборатории попросту не могут се-

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

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

фициентов более высокого порядка. В частности,

ток подобного подхода заключается в том, что экс-

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

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

формирования и положение дисперсионной волны,

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

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

в котором могут проводиться измерения.

спектральной ширины гребенки

[37], а также

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

положение дополнительных спектральных компо-

ровки лишен перечисленных недостатков. Он уни-

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

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

четырехволнового взаимодействия [38, 39].

любого диапазона длин волн, а что самое главное

Пропуская через микрорезонатор излучение от

он достаточно прост в реализации и требует сущест-

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

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

мощность на выходе из микрорезонатора, можно по-

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

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

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

очередь, можно определить положения собственных

мости от требований эксперимента и потенциально

частот (мод) микрорезонатора. Однако на практи-

может достигать 10⁴-10⁵ Гц.

ке подобные измерения трудно осуществить, в осо-

В качестве калибровочного инструмента в

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

рассматриваемом методе предлагается использо-

новная проблема заключается в том, что характер-

вать классический интерферометр Маха - Цендера

ная величина ОСД (D₁/2π) порядка 10¹⁰-10¹² Гц,

(ИМЦ): одновременная регистрация уровней сиг-

а вклад дисперсии групповых скоростей (D₂/2π) в

нала, прошедшего через откалиброванный ранее

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

интерферометр и исследуемый микрорезонатор,

10⁵-10⁸ Гц.

позволяет с высокой точностью определить межмо-

Таким образом, для измерения дисперсионных

довый интервал микрорезонатора. В данной работе

характеристик необходим источник с широким диа-

измерения проводились в диапазоне 1500-1650 нм

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

и использовался волоконный ИМЦ, собранный из

терагерц) и одновременно достаточно узкой мгно-

двух волоконных делителей 50/50 и одномодового

венной шириной линии (на уровне единиц-десятков

волокна smf-28. ИМЦ должен быть предварительно

килогерц). Стабильность и воспроизводимость час-

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

тоты за время проведения измерений также не

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

должны превышать этой величины.

Схема экспериментальной установки для про-

К сожалению, источников, обладающих одновре-

ведения предварительной калибровки ИМЦ пред-

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

ставлена на рис. 3а. В качестве источника излуче-

дняшний день попросту не существует. Узкополос-

ния использовался перестраиваемый лазер Toptica

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

CTL-1550, на вход которого подавался сигнал тре-

ностью, как правило, или имеют очень малый диа-

угольной формы с импульсного генератора (Ген.

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

имп.) для перестройки частоты. Частота треуголь-

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

ного сигнала перестройки была 10 кГц. Оптическая

применения в рассматриваемом эксперименте. Ис-

схема имела два равнозначных плеча, на выходе

точники с возможностью плавной перестройки в

каждого плеча имелся фотодетектор (ФД), кото-

широком диапазоне, в свою очередь, обычно име-

рый передавал сигнал, пропорциональный мощно-

ют низкую долговременную стабильность. Следо-

сти прошедшего оптического сигнала, на осцилло-

ЖЭТФ, вып. 1 (7)

Н. Ю. Дмитриев, А. С. Волошин, Н. М. Кондратьев и др.

ЖЭТФ, том 162, вып. 1 (7), 2022

23.4,-

#$%&$

+,-./01-

!)*

$%&

'#(

"#"

!"#$%&'$()

"#!

*$+%&,

+,-

$%&

Рис. 3. (В цвете онлайн) Калибровка ИМЦ. а) Принципиальная схема установки для калибровки ИМЦ: Лазер перестра-

иваемый источник лазерного излучения (Toptica CTL 1550); Ген. имп. генератор импульсов произвольной формы для

перестройки частоты лазерного источника (Keysight 33500B); ЭОМ электрооптический модулятор (Thorlabs LN65S-FC);

ВЧ ген. высокочастотный генератор для подачи сигнала на ЭОМ; ВР волоконный резонатор на базе делителя 90/10;

ИМЦ исследуемый волоконный ИМЦ; ФД фотодетектор; ОСЦ осциллограф. б) Данные, записанные с экрана

осциллографа (наличие ЭОМ в одном из плеч приводит к генерации боковых резонансов ВР на частоте, равной частоте

сигнал ВЧ генератора), что позволяет перевести временную шкалу осциллографа в частоту и рассчитать период калибру-

емого ИМЦ. в) Аналитическая зависимость эффективного показателя преломления n_eff используемого волокна smf-28 от

длины волны. г) Зависимость периода калибруемого ИМЦ, полученная путем аппроксимации экспериментальных данных

теоретической зависимостью (2)

граф (ОСЦ). В первом плече располагался калибру-

числяется ее период T_MZI . В процессе постобработ-

емый ИМЦ, выходной сигнал которого представлял

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

собой синусоиду. Во втором плече размещался воло-

ли аппроксимированы функцией Лоренца, а сигнал

конный резонатор (ВР) и фазовый электрооптиче-

ИМЦ синусоидой. Повторив подобный экспери-

ский модулятор (ЭОМ), на вход которого с генера-

мент для 13 разных длин волн из интересующего нас

тора (ВЧ ген.) подавался сигнал с частотой, близкой

диапазона, мы получили зависимость периода ИМЦ

к частоте синусоидального сигнала от исследуемо-

от длины волны (синие точки на рис. 3г). Сред-

го ИМЦ. Волоконный резонатор был собран на базе

неквадратичное отклонение результатов измерений

делителя 90/10, и его межмодовый интервал значи-

периода составило 2 · 10^-6 нм. Используя предостав-

тельно превышал ОСД калибруемого ИМЦ.

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

мощью моделирования методом конечных элемен-

Наличие модуляции во втором плече приводило

тов мы получили аналитическую зависимость эф-

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

фективного показателя (n_eff ) преломления волокна

го резонатора на экране осциллографа можно было

от длины волны (рис. 3в).

наблюдать два боковых резонанса, которые отстоят

от основного точно на частоту сигнала модуляции

(рис. 3б ). Зная это расстояние, можно осуществить

однозначный переход от временной шкалы осцилло-

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

графа к шкале частот и получить частотную зависи-

да ИМЦ для разных длин волн были аппроксимиро-

мость функции пропускания ИМЦ, из которой вы-

ваны выражением (2) для периода ИМЦ с использо-

ЖЭТФ, том 162, вып. 1 (7), 2022

Определение дисперсионных характеристик. . .

ванием аналитической зависимости n_eff (λ) [41,42]:

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

так как используемый лазер с отключенной опцией

T_MZI =

(2)

перестройки позволяет выставить начальную длину

n_g(λ₀)L

волны с точностью 10^-5 нм. Поэтому было доста-

где L разность плеч ИМЦ, n_g(λ₀) = n_eff (λ₀) -

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

-λ₀dn_eff/dλ. В результате аппроксимации была по-

ной опцией перестройки и начать сбор данных до

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

включения перестройки. В общем же случае, при

ИМЦ от длины волны T_MZI (λ) (красная линия на

отсутствии описанной возможности лазера, можно

рис. 3г). Таким образом, была проведена калибров-

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

ка ИМЦ. С помощью полученной зависимости мож-

оптическую схему. В этом плече можно осуществить

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

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

дый момент времени. Для этого требуется знать зна-

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

чение частоты источника лишь в одной точке.

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

Следующим этапом эксперимента было непо-

ными спектральными характеристиками, например,

средственное измерение спектра собственных частот

волоконную брэгговскую решетку.

интегрального резонатора. Принципиальная схе-

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

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

ходных данных к частной заключается в следую-

рис. 4а. В качестве лазера накачки использовался

щем: информация о длине волны λ₀ в какой-либо

тот же лазер Toptica CTL 1550, что и в процессе

одной точке исходного массива данных дает нам

предварительной калибровки ИМЦ. Однако теперь

возможность вычислить период ИМЦ в окрестности

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

этой точки, T_MZI(λ₀), используя ранее полученную

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

зависимость T_MZI(λ). Зная величину T_MZI(λ₀), мы

но, от 1510 до 1620 нм. Оптическая схема имела

можем вычислить длину волны нескольких сосед-

два равнозначных плеча, в конце каждого из кото-

них точек в пределах одного периода ИМЦ. Далее

рых располагался фотодетектор (ФД). Сигнал с ФД

необходимо вычислить период ИМЦ уже для новой

в режиме реального времени записывался с помо-

длины волны, λ₁ = λ₀ + T_MZI (λ₀), используя кото-

щью аналогово-цифрового преобразователя (АЦП),

рый мы можем найти значения длин волн для точек,

конкретно в рассматриваемом эксперименте исполь-

соседних с λ₁. Таким образом, мы имеем рекуррент-

зовался АЦП, встроенный в контроллер Toptica. В

ное соотношение

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

λ_i = λ_i-1 + T_MZI(λ_i-1),

ный ИМЦ, а в другом плече фотонный чип c

интегральным микрорезонатором и волноводом из

с помощью которого можно вычислить длину волны

нитрида кремния (Si₃N₄). Ввод и вывод излучения

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

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

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

щью подведения к торцам чипа конических (линзи-

волны (частотной). Для реализации описанного пе-

рованных) волокон (рис. 4б,в), суммарные потери

рехода были написаны соответствующие оригиналь-

при этом составляли около 3 дБ. Исходный массив

ные скрипты в пакете MATLAB. На рис. 4г приведе-

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

ны спектры пропускания микрорезонатора и ИМЦ

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

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

мощности двух каналов от времени. Сигнал одного

Точность определения положения собственных

из каналов, в котором располагается ИМЦ, пред-

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

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

риода используемого ИМЦ и количества точек,

дом. Сигнал второго канала имеет вид плавно меня-

укладывающегося в этот период. Количество точек

ющейся кривой с характерными узкими провалами,

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

соответствующими модам микрорезонатора.

зации АЦП и скоростью перестройки длины волны

Используя полученную на предыдущем этапе за-

лазерного источника. Оптимальная скорость скани-

висимость периода ИМЦ от длины волны, можно

рования должна быть выбрана из следующих со-

осуществить переход от временной зависимости в

ображений: с одной стороны, скорость сканирова-

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

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

тотной (рис. 4г). Для этого необходимо знать точ-

ты дискретизации АЦП хватало для четкой ¾про-

ную длину волны в какой-либо момент времени. В

рисовки¿ синусоиды ИМЦ и пиков резонатора; с

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

другой стороны, слишком низкая скорость скани-

Н. Ю. Дмитриев, А. С. Волошин, Н. М. Кондратьев и др.

ЖЭТФ, том 162, вып. 1 (7), 2022

!"#$%&'#

!""

()*!+

!""

$%&' %()

!-#

!,#

!"#

!+# "## "-#

!-#

!,#

!"#

!+# "## "-#

#$%

!*#5!*

2,#

2-#

3 (4

2,#

2-#

^- ж

!"#

!"!

!*#

!*!

!+#

!+!

!-#

!,#

!"#

!+# "## "-#

!-#

!,#

!"#

!+# "## "-#

$%&' %()

Рис. 4. (В цвете онлайн) Измерение дисперсии интегрального микрорезонатора. а) Принципиальная схема эксперимен-

тальной установки: ИМЦ ранее откалиброванный волоконный ИМЦ;

фотодетектор; АЦП аналого-цифровой

преобразователь для сбора и записи данных. б) Фото общего вида экспериментальной установки. в) Вид в микроскоп

чипа и подводимых в его торец линзованных волокон (процесс юстировки): показаны интегральные резонаторы с меж-

модовым расстоянием 1 ТГц (вверху) и 150 ГГц (внизу). г) Экспериментальные данные для резонатора с межмодовым

расстоянием 150 ГГц после обработки: синяя линия пропускание интегрального чипа с резонатором (спектр мод ре-

зонатора); красная пропускание ИМЦ (на вставке более крупный масштаб). д-з) Экспериментальные результаты для

резонаторов с межмодовым расстоянием 150 ГГц (д, ж) и 1 ТГц (е, з): д и е зависимости отступа по частоте µ-й

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

дисперсии высших порядков в положения мод и аппроксимация квадратичной функцией

рования будет приводить к тому, что частотная

(рис. 4в). Для дальнейшего расчета коэффициен-

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

тов D₁ и D₂ из соотношения (1) необходимо вы-

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

брать точку разложения моду, на которой будет

фектами, приведет к сильным амплитудным шумам

накачиваться будущая гребенка. Для резонатора с

снимаемых данных. В данной работе эксперимен-

межмодовым интервалом 150 ГГц в качестве моды

тально определенная оптимальная скорость скани-

отсчета µ = 0 была выбрана мода вблизи 1567 нм

рования была около 2 нм/c, период использованно-

(f₀ = 191.32 ТГц), а с интервалом 1 ТГц вблизи

го ИМЦ для 1550 нм составлял 103.412 МГц. Та-

1560 нм (f₀ = 192.17 ТГц). После вычитания часто-

кой выбор обусловлен ограничением максимальной

ты выбранной центральной моды f₀ из положения

частоты дискретизации используемого АЦП, кото-

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

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

µ-й моды от центральной моды вида (рис. 4д, е)

чить около 10-12 точек на период синусоиды ИМЦ

(T_MZI (λ)).

∆f(µ) = f_µ - f₀,

а в соответствии с соотношением (1)

2.2. Результаты

D₁

D₂

f_µ - f₀ =

µ+

µ².

2π

4π

На описанной экспериментальной установке с ис-

пользованием изложенной методики нами успеш-

Следовательно, для вычисления коэффициентов D₁

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

и D₂ достаточноаппроксимироватьполученную экс-

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

периментальную зависимость квадратичной функ-

резонаторов из нитрида кремния (Si₃N₄) с меж-

цией. В связи с тем, что D₁ ≫ D₂, удобнее сначала

модовыми интервалами около 150 ГГц и 1 ТГц

найти D₁ путем аппроксимации экспериментальной

ЖЭТФ, том 162, вып. 1 (7), 2022

Определение дисперсионных характеристик. . .

зависимости f_µ-f₀ линейной функцией, а затем вы-

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

числить D₂, найдя квадратичное приближение для

резонатора методом конечных элементов (МКЭ) в

зависимости f_µ - f₀ - (D₁/2π)µ (рис. 4д).

частотной области. В модель была заложена топо-

На графике для резонатора с межмодовым рас-

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

стоянием в 150 ГГц отчетливо видны отклонения

150 ГГц с учетом всех погрешностей геометрических

экспериментальных точек от теоретической кривой,

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

определяемой квадратичной функцией, в районе

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

длин волн 1555, 1595 и 1615 нм. Данные отклоне-

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

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

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

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

Измерение дисперсионных характеристик было про-

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

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

[40,43]. Наличие подобного влияния обусловлено на-

ным направлением вектора напряженности электри-

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

ческого поля (рис. 5а). Именно для этого семейства

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

были рассчитаны положения всех собственных час-

ства. Такие неоднородности уникальны для каждого

тот микрорезонатора в диапазоне 1500-1650 нм. Рас-

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

четы проводились в пакете COMSOL MultiPhysics,

тическом выражении (1).

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

В результате проведенного эксперимента было

ной способностью [45-47]. В результате численного

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

расчета с учетом погрешностей всех параметров для

ладают аномальной дисперсией. Для интегрального

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

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

частот, в которой, согласно теории, она должна рас-

расстоянием 150 ГГц на длине волны 1567 нм

полагаться.

D₁/2π = 143.516 ± 0.005 ГГц,

D₂/2π = 1.72 ± 0.02 МГц,

Таким образом, теоретически был получен до-

а для межмодового расстояния 1 ТГц на длине вол-

верительный интервал для коэффициента D₂/2π =

ны 1560 нм

= 1.11-2.14 МГц (красная область на рис. 5б). Экс-

D₁/2π = 996.80 ± 0.07 ГГц,

периментально измеренное значение коэффициента

D₂/2π = 1.72±0.02 МГц находится в хорошем соот-

D₂/2π = 14.9 ± 0.8 МГц.

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

Погрешность расчета показателей D₁ и D₂ была

оценена с помощью метода складного ножа [44],

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

В качестве дополнительной верификации было

одной экспериментальной точки из расчета. Вели-

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

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

ментальных результатов с уже зарекомендовавшим

межмодовым расстоянием 1 ТГц выше, так как ко-

себя методом измерения дисперсионных характери-

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

стик с использованием источника внешней эталон-

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

ной частотной гребенки, обладающим высокой точ-

с межмодовым расстоянием 150 ГГц (10 мод в пер-

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

вом случае и порядка 70 мод во втором). Предлага-

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

емая методика может быть доработана для расчета

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

дисперсионных коэффициентов более высокого по-

зонаторов использованных в работе [48], где калиб-

рядка.

ровка в процессе измерения проводилась с помощью

источника внешней оптической частотной гребенки,

2.3. Верификация результатов измерения

используемого в качестве эталонного источника. Ко-

дисперсионных характеристик

эффициенты D₂/2π для данной топологии микроре-

Верификация описанного метода была проведе-

зонатора, найденные с помощью предлагаемого ме-

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

тода и приведенные в работе [48], составили соответ-

стоянием 150 ГГц путем сравнения полученных экс-

ственно 1.26 МГц и 1.56 МГц. Отклонение резуль-

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

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

ленных расчетов.

точностью процесса производства.

Н. Ю. Дмитриев, А. С. Волошин, Н. М. Кондратьев и др.

ЖЭТФ, том 162, вып. 1 (7), 2022

;'1(% &#'(!* (+

<3"=821+8($

>%"?8$

0121(% &#'(#&#)%* +3+

7#+82 +#)! 9

Рис. 5. (В цвете онлайн) Сравнение полученных экспериментальных результатов с результатами численного расчета

с помощью метода конечных элементов в частотной области. а) Распределение напряженности электрического по-

ля в поперечном сечении волновода для центральной моды f0 исследуемого пространственного семейства. б) Вклад

дисперсии высших порядков в положения мод. Синие точки

эксперимент, оценка D2/2π = 1.72 МГц. Красная

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

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

D2/2π = 1.11-2.14 МГц

3. ГЕНЕРАЦИЯ СОЛИТОННОЙ

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

ЧАСТОТНОЙ ГРЕБЕНКИ В

ции [57-60].

ИССЛЕДУЕМЫХ МИКРОРЕЗОНАТОРАХ

Излучение лазерного диода в чип заводилось че-

рез торец (рис. 6 б ). С этой целью была изготов-

Солитонная частотная гребенка возбуждалась

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

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

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

ванных микрорезонаторов с D₁/2 = 143.516 ГГц,

прецизионной шестикоординатной подаче. Как из-

D₂/2π = 1.72 МГц, добротностью Q ≈ 2 · 10⁶, соб-

вестно, эффект затягивания является крайне чув-

ственной шириной линии κ₀/2π = 106 МГц и шири-

ствительным к фазе отраженной волны [11], поэто-

ной линии, обусловленной связью, κ_c/2π = 75 МГц

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

(рис. 6 в). В качестве источника накачки использо-

ода в плоскости фотонного чипа и по вертикали к

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

нему требовалось с высокой точностью контролиро-

с резонатором Фабри - Перо, работающий в режиме

вать расстояние между торцами фотонного чипа и

затягивания [2,11]. Работа в этом режиме, обеспечи-

лазерного диода. Для этого были использованы пье-

вающем за счет обратного рэлеевского рассеяния на

зоэлементы, позволяющие позиционировать диод с

поверхностных и объемных неоднородностях микро-

точностью до 20 нм. Чип при этом был термостаб-

резонатора [49, 50] быструю обратную связь между

лилизован и располагался неподвижно на специаль-

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

но изготовленном держателе. К противоположно-

но стабилизировать частоту генерации и значитель-

му торцу фотонного чипа было подведено кониче-

но уменьшить ширину линии [2,30], причем не толь-

ское волокно. Выходной сигнал подавался на опти-

ко для одночастотных, но и для многочастотных ди-

ческий анализатор спектра и детектор, подключен-

одов [51]. Это дает возможность использовать для

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

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

лографу (рис. 6 а). Таким образом, собранный экс-

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

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

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

лучать информацию об оптическом спектре и элек-

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

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

ные диоды [11, 48, 52-54]. Более того, в случае гене-

оценивать степень его когерентности. Как известно,

рации гребенки при помощи затянутого лазерного

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

диода происходит компенсация нежелательных теп-

увеличением тока. Таким образом, возбуждение ин-

ловых эффектов [55,56], которые неизбежны при ге-

тересующей нас моды микрорезонатора и затягива-

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

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

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

ного увеличения тока последнего. При попадании в

ЖЭТФ, том 162, вып. 1 (7), 2022

Определение дисперсионных характеристик. . .

область затягивания сначала наблюдается снижение

ски [32]. Отметим, что мощность линий гребенки (за

выходной мощности и характерное изменение опти-

исключением центральной) не зависит от мощности

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

накачки при фиксированной отстройке.

В случае многочастотного лазерного диода с ре-

Таким образом, дисперсия может быть измере-

зонатором Фабри- Перо наблюдается переход к од-

на как по ширине гребенки (числу линий по уров-

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

ню 3 дБ),

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

√

генерации [48,51,61]. Затем при условии оптималь-

(√ )

2(ω₀ - ω)

ной фазы обратной волны путем дальнейшей под-

N_comb =

2 arcosh

D₂

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

наблюдать гиперпараметрическую генерацию. Как

так и по суммарной мощности боковых линий (без

правило, при увеличении отстройки сначала наблю-

учета центральной линии),

дается ¾шумный¿ режим, затем бризерные состоя-

)

ния, и уже после можно достичь когерентного соли-

√

4η²P_th

тонного состояния, которое характеризуется отсут-

P_lines =

2(ω₀ - ω)D₂ - D₂

ствием низкочастотных составляющих (ниже часто-

ты повторения солитонов) [32, 61]. Подобный спо-

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

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

нове измеренных дисперсионных коэффициентов

ра может быть реализован благодаря наличию эф-

огибающей частотной гребенки и эксперименталь-

фекта затягивания, обеспечивающим дополнитель-

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

ную стабилизацию системы лазер-микрорезонатор.

рис. 6г. Ширина теоретически предсказанной гре-

В случае использования для накачки оптически изо-

бенки практически полностью совпала с шириной

лированного внешнего лазера генерация солитона

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

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

ет о точности измеренных дисперсионных характе-

ставляется возможной и необходимо использовать

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

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

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

модуляцию мощности [57] накачки для компенсации

Однако при этом спектр наблюдаемой гребенки вы-

влияния нежелательных тепловых эффектов.

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

На собранной экспериментальной установке с по-

кой расчетной огибающей, что приводит к тому,

мощью описанной выше методики мы успешно про-

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

наблюдали генерацию солитонной частотной опти-

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

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

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

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

рядом причин.

да с резонатором Фабри - Перо мог варьироваться в

Во-первых, вероятно, что мы наблюдали мно-

диапазоне 0-500 мА. Солитонная гребенка наблюда-

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

лась на токах около 300 мА, что соответствует при-

существует не один, а несколько солитонов одно-

мерно 30 мВт в интегральном волноводе, при этом

временно. Подобное состояние вполне может быть

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

устойчивым и когерентным, т. е. с точки зрения низ-

(P_th) для данного микрорезонатора равен 14 мВт.

кочастотных шумов многосолитонное и односоли-

Как было упомянуто ранее, значение коэффици-

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

ента D₂ напрямую влияет на спектральную ширину

состояние как раз характеризуется изрезанной оги-

генерируемой гребенки [32]. Ее огибающая в случае

бающей оптического спектра, причем степень изре-

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

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

выражением:

ния солитонов внутри микрорезонатора [32].

(

√

)

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

D₂

πµ

D₂

P_µ ∼ 2η²Pth

sech²

(3)

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

2(ω₀ - ω)

ром Фабри - Перо вполне возможно наличие Вернье

где η = κ_c/(κ₀ +κ_c) нормированный коэффициент

эффекта, когда некоторые подавленные продольные

связи [0;1), κ = κ_c + κ₀ ширина линии нагружен-

моды лазерного диода совпадают с соответствующи-

ного микрорезонатора. Поправка к мощности цен-

ми модами резонатора (аналог многочастотной на-

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

качки). В этом случае может происходить усиление

кости, но также может быть рассчитана аналитиче-

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

Н. Ю. Дмитриев, А. С. Волошин, Н. М. Кондратьев и др.

ЖЭТФ, том 162, вып. 1 (7), 2022

*+,-".#/-

%&'

*+,-./ .0

1./ 2

' "

5,00778

'#"

')"

'$"

'!"

*+,-3

4--/3

'%"

'("

!!"#$#

5,00 ""6789:.

!""

!#"

!$"

!%"

!&"

%""

Рис. 6. (В цвете онлайн) Генерация оптической частотной гребенки в режиме затягивания. а) Принципиальная схема экс-

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

диод; КД контроллер лазерного диода (ток инжекции и термостабилизация); БФД быстрый фотодетектор; ЭСА

электрический анализатор спектра; ОСА оптический анализатор спектра. б) Фото лазерного диода (ЛД), подведенного

в торец фотонного чипа (SiN); вид в микроскоп. в) Экспериментально измеренная накачиваемая мода микрорезонатора

в линейном режиме (синяя линия) и аппроксимация функцией Лоренца с учетом связи (красная штриховая линия) и

ее параметры. г) Экспериментально полученная солитонная частотная гребенка (синие линии) и огибающая (3), рассчи-

танная на основе измеренных параметров микрорезонатора (красная линия). Вставка: СВЧ-спектр выходного сигнала.

Отсутствие низкочастотных составляющих свидетельствует о высокой когерентности выходного сигнала

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

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

се его усиления. Аналогичный эффект был недавно

ристик интегральных микрорезонаторов с большим

продемонстрирвон для генерации частотных гребе-

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

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

продемонстрировано на примере измерения ДГС

нормальной ДГС [62].

высокодобротных микрорезонаторов из нитрида

В-третьих, изрезанность огибающей оптической

кремния с межмодовыми расстояниями 1 ТГц и

гребенки может быть вызвана наличием неоднород-

150 ГГц. Точность метода была подтверждена как

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

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

действию мод разных семейств, наличие которых

с результатами, полученными численно с помощью

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

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

но измеренной дисперсионной кривой (см. рис. 5б ).

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

лазерным диодом солитонной гребенки при эф-

фекте затягивания. Параметры огибающей спектра

4. ВЫВОДЫ

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

ской оценкой на основе измеренных дисперсионных

Таким образом, авторами разработан и испытан

параметров микрорезонатора.

оригинальный метод экспериментального опреде-

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

Благодарности. Работа выполнена с ис-

микрорезонаторов и подробно описан способ его

пользованием оборудования ЦКП ВНИИОФИ

реализации. Рассмотренный метод, основанный

(ckp.vniiofi.ru) и ЦКП

¾Визуализации высокого

на применении интерферометра Маха - Цендера

разрешения¿ Сколтеха (https://www.skoltech.ru).

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

Финансирование. Работа выполнена при фи-

чается своей универсальностью с точки зрения

нансовой поддержке Российского научного фонда

диапазона измерений, доступностью и высоким

(проект № 21-72-00132).

уровнем точности. Помимо этого, данный метод

ЖЭТФ, том 162, вып. 1 (7), 2022

Определение дисперсионных характеристик. . .

ЛИТЕРАТУРА

22.

P. Marin-Palomo, J. N. Kemal, M. Karpov et al., Na-

ture 546, 274 (2017).

V. B. Braginsky, M. L. Gorodetsky, and V. S. Ilchen-

ko, Phys. Lett. A 137, 393 (1989).

23.

A. Fülöp, M. Mazur, A. Lorences-Riesgo et al., Na-

ture Comm. 9, 1598 (2018).

V. V. Vassiliev, V. L. Velichansky, V. S. Ilchenko et

24.

M. G. Suh and K. J. Vahala, Science 359, 884 (2018).

al., Opt. Comm. 158, 305 (1998).

25.

P. Trocha, D. Ganin, M. Karpov et al., Science 359,

K. J. Vahala, Nature 424, 839 (2003).

887 (2018).

A. B. Matsko and V. S. Ilchenko, IEEE J. Sel. Top.

26.

Q. F. Yang, B. Shen, H. Wang et al., Science 363,

Quant. Electron. 12, 3 (2006).

965 (2019).

A. B. Matsko and V. S. Ilchenko, IEEE J. Sel. Top.

27.

E. Obrzud, M. Rainer, A. Harutyunyan et al., Nature

Quant. Electron. 12, 15 (2006).

Photonics 13, 31 (2019).

I. S. Grudinin, A. B. Matsko, A. A. Savchenkov et

28.

M.-G. Suh, X. Yi, Y.-H. Lai et al., Nature Photonics

al., Opt. Comm. 265, 33 (2006).

13, 25 (2019).

J. Ward and O. Benson, Laser Photon. Rev. 5, 553

29.

Z. L. Newman, V. Maurice, T. Drake et al., Optica

(2011).

6, 680 (2019).

D. V. Strekalov, C. Marquardt, A. B. Matsko et al.,

30.

W. Liang, V. Ilchenko, A. A. Eliyahu et al., Nature

J. Opt. 18, 123002 (2016).

Comm. 6, 7371 (2015).

G. Lin, A. Coillet, and Y. K. Chembo, Adv. Opt.

31.

D. T. Spencer, T. Drake, T. C. Briles et al., Nature

Photonics 9, 828 (2017).

557, 81 (2018).

32.

T. Herr, V. Brasch, J. D. Jost et al., Nature Photon-

10.

M. Kues, C. Reimer, J. M. Lukens et al., Nature Pho-

ics 8, 145 (2014).

tonics 1, 170 (2019).

33.

T. J. Kippenberg, A. L. Gaeta, M. Lipson, and

11.

N. M. Kondratiev, V. E. Lobanov, A. V. Cherenkov

M. L. Gorodesky, Science 361, eaan8083 (2018).

et al., Opt. Express 25, 28167 (2017).

34.

C. Godey, I. V. Balakireva, A. Coillet, and

12.

P. Del’Haye, A. Schliesser, O. Arcizet et al., Nature

Y. K. Chembo, Phys. Rev. A 89, 063814 (2014).

450, 1214 (2007).

35.

S. Fujii and T. Tanabe, Nanophotonics 9,

1087

13.

T. J. Kippenberg, S. M. Spillane, and K. J. Vahala,

(2020).

Phys. Rev. Lett. 93, 083904 (2004).

36.

C. Xu, J. Ma, Ch. Ke et al., Appl. Phys. Lett. 114,

14.

T. J. Kippenberg, R. Holzwarth, and S. A. Diddams,

091104 (2019).

Science 332, 555 (2011).

37.

V. Brasch, M. Geiselmann, T. Herr et al., Science

15.

A. A. Savchenkov, A. B. Matsko, D. Strekalov et al.,

351, 357 (2016).

Phys. Rev. Lett. 93, 243905 (2004).

38.

N. L. B. Sayson, K. E. Webb, S. Coen et al., Opt.

16.

T. Herr, J. Riemensberger, C. Wang et al., Nature

Lett. 42, 5190 (2017).

Photonics 6, 480 (2012).

39.

N. L. B. Sayson, T. Bi, V. Ng et al., Nature Photonics

17.

A. Pasquazi, M. Peccianti, L. Razzari et al., Phys.

13, 701 (2019).

Rep. 729, 1 (2018).

40.

P. Del’Haye, O. Arcizet, M. L. Gorodetsky et al., Na-

ture Photonics 3(9), 529 (2009).

18.

A. L. Gaeta, M. Lipson, and T. J. Kippenberg, Na-

ture Photonics 13, 158 (2019).

41.

L. C. Schneider, Bloomington (2017).

19.

Y. K. Chembo, Nanophotonics 5, 214 (2016).

42.

Z. Tian, S. S.-H. Yam, J. Barnes et al., IEEE Photon.

Techn. Lett. 20, 626 (2008).

20.

T. Udem, R. Holzwarth, and T. W. Hänsch, Nature

416, 233 (2002).

43.

T. Herr, V. Brasch, J. D. Jost et al., Phys. Rev. Lett.

113, 123901 (2014).

21.

T. Fortier and E. Baumann, Comm. Phys.

2, 153

(2019).

44.

B. Efron and C. Stein, Ann. Stat. 9, 586 (1981).