Численное исследование свойств решений нелинейного уравнения Шредингера при распространении лазерных импульсов в световодах Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 13, № 6, 2008

Численное исследование свойств решений нелинейного уравнения Шредингера при распространении лазерных

>к

импульсов в световодах

В. Э. Витковский, М. П. Федорук Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]

We present numerical simulations of nonlinear Schrodinger (NLS) equation with cubic nonlinearity for the Gauss and ring (m = 1) initial spatial profiles. We determine a threshold power for the self-focusing collapse in a bulk dielectric medium and analyze solution properties that are independent of geometry of the experiment and the initial conditions. We observe oscillatory behaviour for different magnitudes near critical power and phase singularities. To approximate smooth solutions of this problem we construct special boundary conditions and implement a numerical algorithm based on a parallel sweep method for implicit Crank—Nicolson scheme. We then equip this scheme with an adaptive spatial and temporal mesh refinement mechanism that enables the numerical technique to correctly approximate singular solutions of the NLS equation. The calculations were performed using a high performance computer cluster allowing acceleration of 28 times over the sequential algorithm.

Введение

В данной работе численно исследуются свойства решений нелинейного уравнения Шредингера (НУШ) с кубической нелинейностью в цилиндрических координатах для начальных условий (г = 0) типа “чирпованых” гауссова импульса Фо (в дальнейшем обозначается как О) и кольцевого распределения с т =1 Фят (К1). Параметр чирпа в представленных ниже расчетах С = 50, что соответствует фокусному расстоянию около 2280 мкм. Входной радиус лазерного пучка а = 100 мкм (на поверхности образца).

Уравнение Шредингера с кубической нелинейностью описывает эффект Керра, заключающийся в зависимости показателя преломления среды от квадрата модуля электрического поля светового излучения [1]. Таким образом, в физической постановке моделируется нелинейное распространение лазерного импульса в прозрачной диэлектрической среде [2] в приближении медленной огибающей электрического поля Ф:

дФ 1 2 д2 д2

’дг + 2Пал±ф + к№|Ф| Ф = °’ А± = дх2 + ду2■ (1)

*Работа выполнена при поддержке междисциплинарного интеграционного проекта СО РАН № 31. © Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

Для данного уравнения безразмерные переменные и единицы измерения используемых в расчетах величин записываются в следующем виде:

z = r = —, Ld = 2nokor0, ko = Г0 = 1^m

L

D

Г0

Л0

Л0 = G.8^m, n0 = 1.4533, n2 = 2.66 ■ 10

-16

cm

"W

I(r', z') = |Ф(г', z')|2

Л0 |A(r,z)|2,

8n2n0n2r^

P

Л2

Л0

8n2n0n2 I0 =

|A|2dr, Po = 22ПЛ0 = 1.32MW,

8n2n0n2

Л20

8n2n0n2r5

2G.9678

TW

cm2

(2)

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

дФ . т |т|2т Л 1 д / дФ

г— + ArФ + |Ф|2Ф = 0, Ar = - — г—

дz г дг V дг

(3)

с начальными условиями вида

Фо(г) = Ag exp (- (1 +)r ) , Фят(г) = ARmrm exp (- (1 +)r ) , m = 1. (4)

2а2

2а2

Параметр чирпа — С = —(—------------)Д—. Здесь — — расстояние от линзы до точки

АоПо(1 - П2)

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

v/a2T72

численная апертура,

А0 — длина волны лазерного излучения. В работе за фокусное расстояние принимается величина f = — — й.

а

1. Интегральная теория НУШ

Следствия теоремы Таланова предсказывают параболическую зависимость среднеквадратичного радиуса (СКР) светового пучка от пройденного пучком расстояния г; стремление СКР к нулю при критической мощности, величина которой не зависит от маш-табных характеристик пучка (радиуса и фокусного расстояния [3, 4]). Интегралами движения для НУШ являются гамильтониан и мощность:

H

дФ^, г)

дг

гdr; P = |Ф^,г)|2 гdr.

(5)

С

С

2

Ниже приведены следующие из теоремы о вириале выражения для среднеквадратичного радиуса и его фокусного расстояния для гауссова и кольцевого распределений соответственно [5]:

1 — Р?) ка2С

р 5 ШІП р ^

1 + С 2 - рс 1 + С - ре

р ег р ег

а2 I 1 - р

(6)

Р^М ка2С

(пЯ1 )2 _ V рег / ^Ш — ка с

(пшіп) р ,

Р ^ ШІП р

1 + с2 - 1 + с2 -

рЯ1 рЯ1

СГ СГ

Верхняя граница критической мощности определяется из условия равенства нулю гамильтониана [6, 7]:

СО

2/ ^ |2г^г •/ |Vf |2г^г

Р/г = --------оо----------------, Р? = 2, Р*1 = 4. (7)

/ |f |4^г о

Нижняя граница критической мощности для НУШ следует из фундаментального, убывающего на бесконечности стационарного решения НУШ (мода Таунса) [8, 9]:

Д±Я — Я + Я3 = 0, Я'(0) = 0, Я(то) = 0,

(8)

^СГ / г2 \

п(г) ~ л/ ПВ2ех^-’ в2 -0'8, ^ег - 11-67, рСТ _1-86-

2. Численный алгоритм

Численные расчеты проводились по разностной схеме Кранка—Николсон а пользовалась сетка с переменным шагом по г:

0.5, ис-

.Фп - Фп

+ аЬФп + (1 — ^)^Фп + |Фп|2(^Ф п + (1 — ^)Фп) _ 0,

ЬФп _ —

1

Гп

Гп +

к

п+1

Фп+1 — Фп

к,

п+1

кп\ ( Фп Фп— 1

кп

к к

кп кп+1

т "Г"

к

кп+1 _ кп + м, м _ 25000, к1 _ 10

—6

(9)

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

2

2

Граничные условия в центре пучка (слева) формулируются исходя из разностной схемы НУШ (9), построенной в декартовых координатах для Ф(^, X, У) [10]:

. Фс - Фо , 0 (Фі - Фо) , 0,, , (Фі - Фо) , ,,т, ,2, ,т, , п ^ П

г--------------+ 2а-------;----------- + 2(1 — а)----------------------+ |Фо| (аФо + (1 — а)Фо) = 0,

т

к к1

Ф(0,Хі, 0) = Ф(0, 0,Уі), Хі = Уі = кі.

(10)

По мере приближения к фокусу амплитуда интенсивности электрического поля возрастала более чем на семь порядков (при околокритических значениях мощности) по сравнению с начальной величиной. Вследствие этого, с ростом интенсивности, по ме-

то|Ф(г', 0)|2

, здесь |Ф(г', 0)|2

ре распространения пучка, уменьшался шаг по г: т ( 0)|2

заданное наперед значение интенсивности в центре пучка, шаг начинал уменьшаться, если интенсивность становилась выше этого значения, что значительно увеличивало время расчета.

Точность расчетов определялась по сохранению численных интегралов движения. В расчетах гамильтониан и мощность (5), сохранялись с точностью вплоть до 10-4 %, при распространении по г на два фокусных расстояния с общим числом шагов т до ^ 108.

Для вычислений по разностной схеме использовался алгоритм параллельной прогонки (Н.Н. Яненко, А.Н. Коновалов) [11, 12, 13, 14, 15], реализованный по технологии МРІ на языке Фортран-77. Выполненные тесты по ускорению показали 22-кратное ускорение на 36 процессорах для рабочей сетки, состоящей из 2.5 • 105 ячеек [16] (рис. 1).

Для постановки граничных условий справа при г ^ то (в расчетах г^ = 5а г^ = 7а) для гауссова пучка использовалось точное решение линейного уравнения:

и

которое имеет вид

д Ф Л т

г— + Дг Ф = 0, дг

а

Фс(г,г) = Ас

(11)

(12)

Рис. 1. Ускорение 5 параллельной прогонки относительно последовательного алгоритма: ^ори — число процессоров, штриховая линия 5 = ^ори

Здесь Фо — точное решение НУШ для начального распределения, полученное методом

^2Р

функций Грина; Ао = -------- — нормировочный коэффициент; а = (1 — ВС) — 'В,

а

„ гг2 , ,, П„>2 „ ( (1 + 'С>Г2\ ( 'Т2 В(1 + С2)

в = в +(1 — вс) , в = 02-,5 = ехр(------202^)ехр в

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

Для расчетов с начальным кольцевым распределением на правой границе используется асимптотика решения, полученная из формулы Грина методом перевала при

г ^ то, Ая1 = — -----нормировочный коэффициент:

а2

фЯ1(^,г) = Ая1Г^в ехр ^1 —°вС^ 5 (13)

Этот прием позволяет подавить нефизические возмущения, идущие с правой границы, вплоть до достаточно больших величин г. С другой стороны, при малых значениях эволюционной переменной г возникают нефизические возмущения, распространяющиеся с левой границы из-за особенности начального условия при кольцевом распределении (вторая производная равна бесконечности). Это приводит к необходимости начинать расчеты с очень малым шагом по эволюционной переменной г и с а = 0.5 + 0(к1), пока особенность в центре пучка не исчезает.

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

На графиках зависимости центральной интенсивности от г для начальных распределений при околокритической мощности соответственно Ро = 1.895 и Ряд = 2.131 представлены также кривые /^(г) = —2------, где RFWHм(z) — зависимость ширины пучка на

wнм(г)

полувысоте интенсивности от г (рис. 2, а). Для обоих распределений зависимость /^(г) хорошо описывает кривую 10(г) вблизи фокуса и после фокусировки [17, 18, 19]. Также было обнаруженно, что профиль интенсивности решений НУШ аппроксимируется гауссовым распределением для достаточно больших г:

г , , 2Ре* ( г2

/е"(г-г) = Д---ГГ ехр ( — Д2-----

RFWHM(z) \ RFWHM(z)

с мощностью

Р = RFwнм(z *)10(г*) ей = 2 ’ отличной от входной мощности пучка, по аналогии с нормировочным коэффициентом для распределения Гаусса (12). Здесь /0(г*) — интенсивность в центре пучка г = 0 при любом г * вблизи области фокусировки, RFWHм(z*) — ширина на полувысоте интенсивности. Зависимость 10(г) для кольцевого распределения явно обнаруживает появление второго максимума вслед за главным пиком интенсивности, в то время как для гауссова распределения образование второго максимума только начинается и слабо выражено.

а б

с о г / л 2 • 1.519

Рис. 2. Зависимость интенсивности в центре пучка от г и /е££(г) = —2-—— для гауссова

RF wнм(г)

распределения Р = 1.895 и кольцевого распределения Р = 2.131 (а); кривые зависимости мощности коллапсирующей части РС01 и коэффициента Ре££ от входной мощности Р (б); распределение интенсивности по радиусу Р = 1.8 (в); зависимость Й^НМ(г) Р = 1.8 (г): все зависимости приведены для гауссова и кольцевого распределений; интенсивность I — в единицах /о, г и г — в микрометрах

в

г

Зависимости Ре££ для гауссова и кольцевого распределений от входной мощности представлены на рис. 2, б. Видно, что при стремлении входной мощности к критической Ре££ 1.519 для обоих распределений. На этом графике также представлены кривые

зависимости мощности, содержащейся в коллапсирующей части профиля РС01 для обоих начальных распределений при г = ^™нм. Здесь имеется в виду мощность между центром пучка и ближайшим минимумом, гямй — фокусное расстояние для среднеквадратичного радиуса, ^™нм — фокусное расстояние для ширины на полувысоте интенсивности. В расчетах наблюдалось, что величина гямй для различных Р совпадает с теоретическим значением (6) и отличается от ^™нм. Было замечено, что мощность РСо1 в этой области при стремлении входной мощности к критической стремится к величине

нижней границы коллапса Pcr ^ 1.86. Каждому начальному распределению (гауссово-му или кольцевому) соответствует собственная кривая Peff(P) и Pcoi(P). Зависимости Pcoi(P) и Peff(P) от входной мощности были также рассчитаны для параметров чир-па C = 100 и C = 0. Расчеты обнаружили независимость от фокусного расстояния значений критической мощности, кривых Peff(P) и Pcol(P).

На рис. 2, в показано распределение интенсивности в коллапсирующей части профиля для гауссова и кольцевого распределений. Зависимости Rfwhm(z) для этих распределений (рис. 2, г) имеют явные отличия в соответствии с зависимостью I0(z) (рис. 2, а). На кривой Rfwhm(z), соответствующей кольцевому распределению, проявляется объединение двух пиков из начального профиля в один в виде скачка на расстоянии примерно 2200 мкм.

На рис. 3 показана зависимость максимальной интенсивности (в центре пучка) от мощности входного пучка для гауссова и кольцевого распределений. Здесь обнаруживается одинаковая зависимость от мощности для гауссова и кольцевого распределений. У интенсивности при z = zRMS для кольцевого распределения возникают осцилляции по мощности, схожие по характеру с осцилляциями по мощности у угла наклона Rfwhm(z) после прохождения фокуса. Для кольцевого распределения интенсивность в центре при zRMS приблизительно на два порядка меньше интенсивности при z = zFWHM для мощностей, близких к критической величине. Интересно заметить, что при предкритических мощностях интенсивность в локальном фокусе zFWHM кольцевого распределения на порядок выше интенсивности гауссова распределения.

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

На рис. 4, а, б представлены развертки профилей по мощности в фокусе zFWHM для гауссова и кольцевого распределений соответственно. На рисунках хорошо заметно появление областей более низкой интенсивности по сравнению с окрестностью, вблизи критической мощности. На рисунках с развертками профилей по z для гауссова профиля (рис. 4, в) и кольца (рис. 4, г) в фокальном объеме также возникают световые пустоты. Графики приведены для предкритических мощностей Pq = 1.895 и Pri = 2.131.

Рис. 3. Зависимость интенсивности при г = 0, г = и г = для гауссова и кольцевого

распределений от мощности Р

-10-

15 2000 2500 3000 z

в

Р

1.0Е+03 0 5 1 1 5 2

1.0Е+01..............................................\

1.0Е-01 3.0Е-02 1.0Е-02 3.0Е-03 1.0Е-03 1.0Е-04

-15:

2000

L2500~

г

3000

Рис. 4. Зависимость распределения интенсивности по г от мощности Р при г = ^^нм для гауссова распределения (а); зависимость распределения интенсивности по г от мощности Р при г = %™нм для кольцевого распределения (б); эволюция гауссова профиля Р = 1.895 (в); эволюция кольцевого профиля Р = 2.131 (г)

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

Выводы

В соответствии с результатами работы расчетный профиль интенсивности для начальных распределений хорошо описывается гауссовым профилем с эффективной мощностью Peff для достаточно больших величин z: Ieff(r, z) = —2 ef . exp ( —

rfwhm (z) V rf whm (z)

Рис. 5. Эволюция фазы амплитуды электрического поля: а — гаусс, Р = 1.895; б — кольцо, Р = 2.131

а

Зависимости Peff и Pco\ стремятся к величинам 1.52 и 1.86 соответственно при стремлении входной мощности пучка к критической для каждого из распределений. Каждому начальному распределению соответствуют собственные кривые Peff(P) и Pcoi(P), не зависящие от входного радиуса пучка и фокусного расстояния. Определены критические мощности локального коллапса для распределений PcG = 1.900 и PcR1 = 2.132, при которых центральная часть профиля коллапсирует и интенсивность в ней стремится к бесконечности, а мощность остается постоянной: Pcol = 1.86; среднеквадратичный радиус соответствует теоретическому значению, отличному от нуля. Фокус для среднеквадратичного радиуса zRMS не совпадает с фокусом для ширины на полувысоте интенсивности z = zFWHM. Обнаружены также спиралевидные фазовые особенности на графиках эволюции фазы в координатах z и r.

Выражаем благодарность В.И. Паасонену и Д.Л. Чубарову за многочисленные полезные обсуждения.

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

[1] Уизем Дж. Линейные и нелинейные волны. М.: Мир, 1974. 622 с.

[2] Агравал Г. Нелинейная волоконная оптика. М.: Мир, 1996. 324 с.

[3] Колоколов И.В., Кузнецов Е.А., Мильштейн А.И. и др. Задачи по математическим методам физики. М.: Эдиториал УРСС, 2002. 288 с.

[4] Talanov V. Focusing of light in cubic media // JETP Lett. 1970. Vol. 11. P. 199-201.

[5] Turitsyn S.K., Mezentsev V.K., Fedoruk M.P. et al. Sub-critical regime of

femtosecond inscription // Optics Express. 2007. Vol. 22, N 15. P. 14750.

[6] Johannisson P., Anderson D., Lisak M. et al. Nonlinear Bessel beams //

http://eprintweb.org physics/0305085 (May 2003)

[7] Fibich G., Gaeta A. L. Critical power for self-focusing in bulk media and in hollow waveguides // Opt. Lett. 2000. Vol. 25. P. 335-337.

[8] Chiao R., Garmire E., Townes C. Self-trapping of optical beams // Phys. Rev. Lett. 1964. Vol. 13. P. 479-482.

[9] Weinstein M. Nonlinear Schrodinger equations and sharp interpolationestimates // Commun. Math. Phys. 1983. Vol. 87. P. 567-576.

[10] Паасонен В.И. Граничные условия повышенной точности в полюсах координатных систем // Вычисл. технологии. 2000. Т. 5, № 1. С. 93-105.

[11] Яненко Н.Н., Коновалов А.Н., Бугров А.Н., Шустов Г.В. Об организации параллельных вычислений и “распараллеливание” прогонки // Числ. методы механики сплош. среды. 1978. T. 9, № 7. С. 139-146.

[12] Паасонен В.И. Параллельный алгоритм для компактных схем в неоднородных областях// Вычисл. технологии. 2003. Т. 8, № 3. С. 98-106.

[13] Паасонен В.И. Сходимость параллельного алгоритма для компактных схем в неоднородных областях // Вычисл. технологии. 2005. Т. 10, № 5. С. 81-89.

[14] Воеводин А.Ф. Метод прогонки для разностных уравнений, определенных на комплексе // Журн. вычисл. математики и мат. физики. 1973. T. 13, № 2. С. 494-497.

[15] Кудряшова Т.А., Поляков С.В. О некоторых методах решения краевых задач на многопроцессорных вычислительных системах // Тр. IV междунар. конф. по математическому моделированию. М.: “СТАНКИН”, 2001. Т. 2. С. 134-145.

[16] Витковский В.Э., Федорук М.П. Вычислительная производительность параллельного алгоритма прогонки на кластерных суперкомпьютерах с распределенной памятью // Вычисл. методы и программирование. 2008. Т. 9, № 1. С. 305-310 (http://num-meth.srcc.msu.ru/).

[17] Akrivis G.D., Dougalis V.A., Karakashian O.A., Mckinney W.R. Numerical

approximation of blow-up of radially symmetric solutions of the nonlinear Schroodinger equation// SIAM J. Scientific. 2003. Vol. 25, N 1. P. 186-212.

[18] Fibich G. Adiabatic law for self-focusing of optical beams // Opt. Lett. 1996. Vol. 21.

P. 1735-1737.

[19] Le Mesurier B.J., Christiansen P.L., Gaididei Y.B., Rasmussen J.J. Beam

stabilization in the 2D nonlinear Schrodinger equation with attractive potential by beam

splitting and radiation // Phys. Rev. E. 2004. Vol. 70. P. 046614.

Поступила в редакцию 17 апреля 2008 г., в переработанном виде — 28 июля 2008 г.