Введение

Среди всевозможных режимов распространения лазерного излучения в нелинейной среде, распространение мощного фемтосекундного импульса представляет в настоящий момент особый научный интерес. Области применения таких знаний: дистанционное зондирование, микрофотоника, дистанционное управлении электрическим разрядом [1]. Впервые эксперименты по дальнему распространению фемтосекундных лазерных импульсов были выполнены в середине 1990-ых [2-4]. В этих экспериментах использовался инфракрасный лазер с продолжительностью импульса около 100фс  и мощностью превышающей , т.е. мощностью, достаточной для самофокусировки импульса [5]. В экспериментах наблюдался распад начального пучка на узкие нити длиной несколько метров. Количество возникающих нитей зависело от мощности импульса. В каждой из них была сконцентрирована доля мощности импульса.

В результате самофокусировки растет интенсивность импульса и уменьшается его ширина, но «схлопывания» не происходит из-за дефокусирующего воздействия электронной плазмы, созданной многофотонной ионизацией молекул воздуха. В результате, максимальная интенсивность в филаменте не превышает 10¹⁴Ватт/см² для инфракрасных импульсов. В зоне максимальной интенсивности регистрируется движущийся вдоль оси распространения импульса фокус, след которого принято называть филаментом (от англ. filament – нить), а процесс образования таких структур - филаментация. Возникновение нескольких фокусов при наличии достаточной мощности импульса обычно объясняется мелкомасштабной самофокусировкой шумовых возмущений начального профиля импульса [6] либо отсутствием его аксиальной симметрии [7].

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

В настоящее время четырьмя институтами Франции и Германии выполняется проект «Teramobile» по экспериментальному и численному исследованию распространения мощных фемтосекундных импульсов. В результате этих экспериментов по распространению терраватных импульсов наблюдалось несколько десятков филаментов, которые упорядочивались в группы (кластеры) протяженностью более десяти метров.

Существует значительное количество работ различных авторов, посвященных математическому моделированию процесса филаментации [7;8;14;16]. Рассматриваемая для данного явления система уравнений для медленно меняющейся амплитуды светового поля является нестационарной 3-х мерной задачей. Для того, чтобы иметь возможность сравнить экспериментальные данные с расчетами при условии учета мелкомасштабных возмущений, требуется порядка ~10¹⁶ счетных ячеек. Такое большое количество делает счет уравнений при больших расстояниях слишком долгим. Для качественного исследования образования филамент и их упорядочивания в кластеры применяется упрощенная физическая модель, которая в совокупности с использованием технологий параллельных вычислений позволяет в обозримое время провести счет задачи.

1.    Физическая модель задачи.

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

,  (1.a)

где                                                         (1.б)

и модель Друда [19] для локальной плотности плазмы :

                                               ,                        (1.в)

где λ₀=800нм  - длина волны;  - центральное волновое число; n₂= 3.1·10^-19 см²/Ватт  - индекс преломления эффекта Керра; k^″=0.2фс²/см  - коэффициент дисперсии групповой скорости; τ_K=70фс  - время релаксации; ρ_c=1.8·10²¹см^-3  - критическая плазменная плотность; β^(K)=4.25·10^-98см¹³/Ватт⁷  - коэффициент многофотонного поглощения; - число фотонов, которое требуется для ионизации; σ=5.44·10^-20см²  - коэффициент каскадной ионизации и плазменного поглощения в поперечном сечении для обратного тормозного излучения; σ_K=2.88·10^-99см¹⁶/с·Ватт⁸  - коэффициент многофотонной ионизации; U_i=12.1e Вольт  - промежуточный потенциал ионизации молекул кислорода; ρ_nt=5.4·10¹⁸см^-3  - эффективная плотность нейтральных молекул, равная 20% от стандартной ρ_at=2.7·10¹⁹см^-3 . Рекомбинация электронов в кинетическом уравнении (1.в) не учитывается из-за малой длины рассматриваемых импульсов. Параметр  для моделирования фемтосекундных импульсов полагается равным 0.5, вовлекая в нелинейный член слагаемое, отвечающее за поляризацию молекул.

2.    Усреднение модели.

В целях сокращения времени счета для качественного исследования распространения мощных импульсов в реальных лазерных системах можно применить алгоритм упрощения данной модели, реализованный в [8]. Во-первых, пренебрежимо мало действие обратного тормозного излучения при продолжительностях импульса менее пикосекунды, поэтому положим  (в кинетическом уравнении для локальной плотности плазмы также не учитывается рекомбинация электронов). Также не будем учитывать изменение второй производной по времени , т.к. комплексная амплитуда световой волны медленно меняется со временем. Во-вторых, в предположении, что в области высокой интенсивности излучения световые поля имеют Гауссов профиль по временной координате и могут быть представлены в виде произведения, будем искать решение системы в виде:

,

где  - временной слой, где формируется максимальный по интенсивности пик с полушириной , сохраняющейся все время распространения [8].

Применяя сделанные выше предположения, сначала найдем  из уравнения (1.в):

,         (2)

где . Далее в уравнение 1(a) подставим найденное  и уравнение (1.б), затем умножим получившееся уравнение на  и проинтегрируем по всему временному домену, т.е. . После упрощений получим:

,       (3)

где

, ,

.

В работе [18] было показано влияние комбинационного рассеяния на нелинейный отклик в воздухе при . Также положим значение  ( – длительность импульса) для изучения процесса филаментации при наличии многофотонной ионизации, которая укорачивает длительность импульса. Чтобы обезразмерить уравнение (3), вводим

переменные:

,, ,

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

,  (4)

где параметр  примет значение  при выбранных выше параметрах  и .

Для численного моделирования уравнений данного типа (НУШ) применяются разные методы, призванные уменьшить требуемое для счета процессорное время и сохранить требуемую точность счета.

Спектральный метод позволяет эффективно использовать алгоритм быстрого преобразования Фурье, но реализуется только на равномерных сетках. Преимущество спектрального подхода бесспорно только для простых решений типа движущегося солитона.

Расщепление по процессам - решение уравнения (4) на шаге  заменяется последовательным решением двумерных задач дифракции в линейной среде сначала в плоскостях, параллельных XZ, затем в плоскостях, параллельных YZ, после чего вычисляется нелинейный набег фазы светового поля на том же интервале  в отсутствие дифракции. Для схем расщепления по процессам существенна проблема консервативности, которая требует отдельного обсуждения.

В данной работе использован прямой метод решения уравнения (4), т.е. построена неявная консервативная разностная схема и после применен метод Ньютона.

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

,    (5)

где  – начальная ширина импульса и - начальная амплитуда (безразмерные величины). Степень N обычно берут равной 1 (гауссов импульс) или 2 (супергаусс).

3.    Численая схема. Оценка размерности моделируемой задачи

Решаемая задача представляет собой нелинейное уравнение параболического типа. Уравнение можно считать двумерным (относительно переменных x, y) и нестационарным относительно переменной z, которая в данном случае принимает роль времени. Для решения задачи составим для квадратной области в координатах  равномерную (в общем случае) сетку с шагом  с общим количеством ячеек . Составим неявную разностную схему (значение функции  будем брать в середине j-ой ячейки составленной сетки):

=0,                    ,                                                                             (5)

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

и гамильтониана

.

Нелинейная задача решается методом Ньютона. Условие окончания итерационного процесса в методе Ньютона выбираются следующие:

,

где  - i-ая компонента вектора решения на текущей (s+1)-ой итерации по нелинейности, ε₁ и ε₂  - числа, определяющие необходимую точность сходимости метода.

При решении задачи выбирается область, превышающая ширину моделируемого импульса в 5-6 раз, на границе ставятся условия второго рода . В качестве начальных условий выбирается профиль:

,    (7)

где  – начальная ширина импульса; - начальная амплитуда. Степень N обычно берут равной 1 (гауссов импульс) или 2 (супергаусс).

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

Чтобы оценить необходимые размеры расчетной сетки, воспользуемся следующими условиями и результатами экспериментов [Teramobile]:

·        для изучения процесса филаментации используется мощный (от 100GW до 3TW) лазерный импульс, начальная ширина которого может достигать 2 см;

·        ширина образующейся филаменты, как мелкомасштабного эффекта, постоянна и достигает 100-150 мкм, а интенсивность в фокусе филаменты может превышать интенсивность фона в десятки раз;

·        наблюдаемая длина распространения импульса составляет от десятков до сотен метров.

Таким образом, длина ребра квадрата расчетной сетки может достигать 10 см, чтобы избежать влияния границ. Длина же ребра одной ячейки должна быть не более 10-15 мкм для удовлетворительной передачи характера мелкомасштабных эффектов. Таким образом, количество ячеек в расчетной области  может достигать. Шаг по z выбирается ~, и уменьшается обратно пропорционально максимальной интенсивности . Тем самым обеспечена не только устойчивость, но и точность. Оценка размерности задачи показывает, что для моделирования реального эксперимента по распространению мощного фемточсекундного лазерного импульса на длинные дистанции в атмосфере необходимо использовать методы параллельного программирования.

4.    Результаты вычислений для одиночного Гауссого импульса

Далее перейдем к решению уравнения (4) с учетом всех описанных эффектов. Интерпретация результатов зачастую будет приводиться также в безразмерных величинах.

На рис. 1 показан характер изменения мощности при использовании начального Гауссова профиля с высокой, близкой к пороговой (или максимальной) интенсивностью и мощностью, не на много превышающей критическую ().

Рис.1 Изменение мощности

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

На графиках приведены зависимости по координате распространения  значений мощности и максимальной амплитуды  для Гауссовых профилей с одинаковой начальной амплитудой  и разными полуширинами . Хорошо видно (рис.2-б), что на скорость возрастания пиковой интенсивности влияет не только мощность импульса (при мощности меньше пороговой  имеет место дефокусировка), но и его ширина.

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

Значение гамильтониана  к моменту установления постоянного уровня мощности (рис.2–а) меняет знак на положительный, что, как известно [13], является необходимым условием расходимости пучка в среднем. Однако факт расходимости пучка в среднем не противоречит тому, что некоторая его внутренняя часть может сфокусироваться. Проанализировать дальнейшее развитие профиля можно с помощью теории Беспалова-Таланова. Для этого рассмотрим спектральный состав распределения  на достаточно большом расстоянии  (рис. 3).

Рис. 3. Спектральный состав распределения

Напомним [6], что инкремент неустойчивости играет роль, только если в спектре импульса присутствуют гармоники, меньшие по величине среднего значения распределения , которое для рассматриваемого  почти равно нулю. Поэтому с точки зрения неустойчивости Беспалова-Таланова дальнейшее развитие неустойчивых возмущений не возможно.   

Чтобы объяснить ступенчатый характер падения мощности, рассмотрим подробнее этапы самофокусировки Гауссова импульса, описанные в [9]. Как было сказано ранее, для самофокусировки Гауссова импульса необходима концентрация малых частот спектра. Если ширина пучка велика, то также велика доля низкочастотных компонент и происходит самофокусировка без начального уширения импульса в отличие от того, как это бывает с менее широкими импульсами. В итоге образуется довольно крутой «столб» интенсивности с наличием вокруг фона из дефокусирующейся части начального импульса (рис.4-а). Но при достижении пиковым значением интенсивности порогового значения начинается падение интенсивности как за счет потери энергии на ионизацию, так и дефокусирующего воздействия электронной плазмы, возникшей в результате многофотонной ионизации. Первые участки, на которые оказывается воздействие, являются зонами больших градиентов. В рассматриваемом случае интенсивность падает быстрее на краю «столба» (рис.4-б).

При дальнейшем падении интенсивности отток мощности от центра встречается с областями подпитки и образуется квазистационарная кольцевая структура (рис.4-в). На рис.5 (который соответствует рис.4-б) показаны потоки мощности, определяемые градиентом фазы . На рисунке можно выделить 3 структуры: центральный приосевой профиль и 2 кольцевых структуры. Внутри этих колец поток мощности направлен к центру. Эти структуры можно отделить друг от друга областями с минимальной интенсивностью, где поток мощности направлен «от центра» (заштрихованы на рисунке).

Рис. 5 Градиент фазы распределения комплекснозначной величины .

При концентрации определенной мощности в кольцевой структуре после падения всего «столба» создаются условия для повторного роста интенсивности (рис.4-г). Такими условиями в рассмотренном случае являются форма распределения (рис.4-г) - Гауссов профиль с мощностью, большей критической, и достаточной для развития самофокусировки полушириной. Таким образом, происходит первоначальный рост интенсивности, образование кольцевой структуры и после падения интенсивности формирование такого профиля, который вновь самофокусируется. В итоге получается ступенчатое падение мощности импульса (рис.2-а).

Характер изменения мощности на рис.1 можно объяснить отсутствием вокруг импульса «фона» низкочастотных компонент спектра. Поэтому моделирование распространения импульсов необходимо начинать с момента входа лазерного импульса в нелинейную среду, а не со стадии развитой самофокусировки. Также стоит отметить, что под установлением постоянного уровня мощности (рис.2–а)  понимается лишь сохранение ее до какого то , потому что, как обсуждалось ранее, и в случае устойчивых возмущений возможен процесс филаментации на дальних расстояниях.

На рис.6 одновременно представлены зависимости падения мощности импульса (Power) и возрастания его пиковой интенсивности (Intensity) на протяжении распространения импульса (z). Можно заметить, что уровень, при котором начинается поглощение, равен . Известна оценка [17], что ионизация воздуха становится существенной при интенсивности 4.5·10¹³Ватт/см² , что и соответствует найденному  безразмерному пороговому значению. По мере уменьшения мощности импульса уменьшается величина теряемой мощности и длина, на которой наблюдается поглощение.

­­­

Рис.6 Сопоставление интенсивности и мощности

Напишем уравнение для изменения мощности по . Для этого умножим уравнение (4) на  и сложим с комплексно-сопряженным уравнением для (4), умноженным на . Проинтегрировав по всей поперечной плоскости , получим

; ;

.

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

,

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

Характер ступенчатого развития мощности импульса также зависит от его первоначальной ширины. На рис. 7 приведены графики изменения мощности для импульсов одинаковой мощности. 

Рис.7. Изменение мощности при различных начальных ширинах импульсов ()

5.    Результаты вычислений для двух Гауссовых импульсов

В работе [9] нами уже было рассмотрено взаимное влияние друг на друга двух Гауссовых импульсов. При наличии многофотонной ионизации влияние становится более заметным. Работы по моделированию влияния двух импульсов проводились группой исследователей под руководством В.П. Кандидова и описаны в [10], где приведены условия на мощность импульсов и расстояние между ними, при которых взаимодействие возможно. Но, как и во многих работах, не уделено должное внимание к зависимости соотношения между ширинами импульсов и их интенсивностями – при одинаковых значениях мощности эти соотношения могут широко варьироваться. Кроме того, наглядно не показан характер взаимодействия областей высокой интенсивности импульсов и не рассмотрен вопрос о точке () начала взаимодействия.

Приведем два показательных примера расчета двух Гауссовых профилей, демонстрирующих, как наличие диссипации влияет на характер развития областей высокой интенсивности. Были рассмотрены два равных Гауссовых возмущений с полушириной  и амплитудой интенсивности  на расстоянии . Распространение импульса моделировалось на расстояние .

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

Рис.9 Переток мощности в широком пучке.

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

6.    Результаты расчета мощного лазерного импульса

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

Мощность реального лазерного импульса достигает нескольких терраватт, ширина импульса может достигать нескольких сантиметров. Профиль импульса реальных лазерных систем обычно задают в форме функции супер-гаусса, т.е. показатель степени  в формуле (5). Особенностью реального лазерного пучка также является наличие шума во входном профиле. Для воспроизведения случайного поля флуктуаций амплитуды использовался спектральный метод [15]. Этот метод основывается на суммировании Фурье-гармоник пространственного спектра. Алгоритм задания спектральным методом случайной реализации комплексного поля флуктуаций входного сигнала состоит в следующем. Сначала получаем комплексную случайную величину для каждой точки спектрального пространства , где  и - независимые равномерно распределенные на отрезке [-1,1] случайные величины. Затем строится функция

, где

- амплитуда шума (например, 0.1),  - радиус корреляции (например, 1/0.1мм). Для получения амплитудного шума выполняется преобразование Фурье:

Окончательно, накладываем функцию шума на входной профиль по следующему закону

.

Развитие мелкомасштабной самофокусировки [6], приводящей к распаду симметрии профиля, классически объясняется наличием амплитудных флуктуаций начального распределения. Вследствие процесса мелкомасштабной самофокусировки и достижения амплитудами отдельных флуктуаций порога многофотонной ионизации в мощном лазерном импульсе начинает развиваться множественная филаментация. При наблюдении этого процесса вырисовывается картина множественных связей отдельных филамент между собой, вследствие чего происходят поглощение и зарождение новых филамент. В процессе такого хаотического развития в лазерных импульсах можно выделить области, где процесс филаментации продолжается непрерывно за счет поглощения энергии соседних участков. Такие области называют кластерами [1] и именно их наличие объясняет характер экспериментальных данных (рис.10).

 Чтобы оценить необходимые размеры расчетной сетки, воспользуемся следующими условиями и результатами экспериментов (проект Teramobile [8]):

·        для изучения процесса филаментации используется мощный (от 100GW до 3TW) лазерный импульс, начальная ширина которого может достигать 2 см;

·        ширина образующейся филаменты, как мелкомасштабного эффекта, постоянна и достигает 100-150 мкм, а интенсивность в фокусе филаменты может превышать интенсивность фона в десятки раз;

·        наблюдаемая длина распространения импульса составляет от десятков до сотен метров.

Таким образом, длина ребра квадрата расчетной сетки может достигать 10 см, чтобы избежать существенного влияния границ. Длина же ребра одной ячейки должна быть не более 10-15 мкм для удовлетворительной передачи характера мелкомасштабных эффектов. Таким образом, количество ячеек в расчетной области  может достигать . Шаг по  выбирается ~, и уменьшается обратно пропорционально максимальной интенсивности . Тем самым обеспечена не только устойчивость, но и точность. Ввиду большой размерности решаемой системы линейных уравнений и необходимости численного моделирования процесса распространения фемтосекундного импульса на длинные дистанции использовался параллельный вычислительный комплекс.

Рис.10 Эксперимент по распространению лазерного луча

В качестве примера расчета приводим расчет для среднего по мощности () импульса с полушириной ˜ω₀= 3мм . Расчетная область представляет собой квадрат с длиной стороны 2см . На начальном профиле задан Гауссов шум по описанному выше алгоритму. Расчет проводился на сетке  ячеек.

Рис.12 Линии уровня распространения импульса

с мощностью  на расстояние около 9м .

Результаты развития мощности (рис. 13-а) и интенсивности (рис. 11,12,13-б) импульса совпадают с полученными результатами проекта Teramobile [8], более того сравнение временных характеристик по скорости счета созданного для данной работы параллельного комплекса также дает хорошие результаты, не смотря на коренное различие с [8] в подходах к численному решению задачи. Подробное обсуждение используемых методов аппроксимации и характеристик параллельного алгоритма приводится в [20]. Как правило, для большей кучности филамент, на практике применяется эллиптическое распределение начального профиля.

Заключение

Целью настоящей работы было описание свойств процесса филаментации и разработка эффективных алгоритмов для расчета этого явления в реальных экспериментальных лазерных системах. Были использованы симметричные разностные схемы, обладающие свойством консервативности, созданы распараллеленные программы для расчета таких задач.

Показан ступенчатый характер развития нелинейности, который определяется пороговым значением интенсивности 5·10¹³Ватт/см² , при котором ионизация воздуха становится существенной. Характер процесса определяется соотношением между процессом многофотонной ионизации и мелкомасштабной самофокусировки. До тех пор, пока интенсивность не достигает пороговых значений, световой пучок не теряет энергию и происходит самофокусировка. При достижении пороговых значений интенсивности начинается потеря энергии на ионизацию и дефокусирующее воздействие электронной плазмы, что приводит к потере интенсивности. Когда интенсивность падает ниже порогового значения, вновь может начаться самофокусировка. Таким образом, процесс носит ступенчатый характер (рис. 10). Процесс прекращается, когда спектральный состав становится таким, что дальнейшее развитие неустойчивости в линейном приближении невозможно, т.к. выполняется условие устойчивости согласно критерию Беспалова-Таланова. Однако, не доказана возможность его возобновления на протяжении большого расстояния.

11. Aztec. A Massively Parallel Iterative Solver Library for Solving Sparse Linear Systems (Инетернет-ресурс: http://www.cs.sandia.gov/CRF/aztec1.html)

12. Дацюк В.Н., Букатов А.А., Жегуло А.И.. Многопроцессорные системы и параллельное программирование (Ростов: Ростовский государственный университет, http://rsusu1.rnd.runnet.ru/koi8/index.html)

13. Ландау Л.Д., Лившиц Е.М. Электродинамика сплошных сред (М.: Наука, 1959).

14. Карамзин Ю.Н., Сухоруков А.П., Трофимов В.А. Математическое моделирование в нелинейной оптике (М.: Издательство Московского университета, 1989).

15. Миркин Л.И., Рабинович М.А., Ярославский Л.П. ЖВМ и МФ, 5, 1353 (1972)

16. Волков В.М. Дифференциальные уравнения, 29, 7 (1993)

17. Kasparian J., Sauerbrey R., Chin S.L. Appl. Phys. B, 71, 877 (2000)

18. Nibbering E.T.J., Grillon G., Franco M.A. at al. J.Opt.SocAm.B, 14, 650 (1997)

19. Feit M.D., Fleck J.A. Appl. Phys. Lett., 24, 169 (1974)

20. Балашов А.Д., Пергамент А.Х. Математическое моделирование, 18, 4, (2006)