Glava1.tex

\newpage
\clearpage
\section{Результаты работы}

\hspace{0.5cm} Результаты исследовательской практики, представленные в настоящей работе, были получены с 10.2014 по 06.2016 в лаборатории вычислительной аэродинамики ИТПМ СО РАН и лаборатории неравновесных течений НГУ.

Моделирование проводилось с использованием программных комплексов $Smile$ (Fortran77) и $Smile$++ (C++). Первые результаты, в том числе имплементация аналитического решения \cite{ShnKinetic} затухания интенсивности в возмущенном газе, были получены на $Smile$. Было проведено их сравнение с результатами более ранних работ M. Shneider, S. Gimelshein, T. Lilly, J. Graul и др. Дальнейшее развитие численной модели и основные результаты, представленные в данной главе, получены на $Smile++$ - более современной, гибкой и удобной для разработке платформе, обладающей также рядом значительных улучшений в производительности и скорости вычислений по сравнению с ранними вычислительными кодами.

В ходе работы были получены следующие результаты:

\begin{itemize}
\item Разработан новой программный модуль, описывающий воздействие оптической решетки как внешней силы на газ поляризуемых частиц, в том числе с учетом обратного взаимодействия газа на излучения на основе аналитических моделей.
\item Создана и верифицирована модель тонких слоев \cite{Born}, описывающая прохождение излучения через газ как оптическую среду, решая полную систему уравнений Максвелла для электромагнитного поля.
\item Исследовано обратное влияние газа на излучение для различных параметров взаимодействия - размера области взаимодействия, свойств газа и излучения, в том числе на вычислительных кластерах.
\item Проведено сравнение результатов, выполнен анализ полученных данных.
\end{itemize}

Предшествующие исследования, на которые опирается данная работа, были выполнены (за некоторыми исключениями) в предположении, что обратное влияние газа на излучение слишком мало или отсутствует. Это предположение, как мы увидим позднее, обоснованно, если область взаимодействия, интенсивность излучения, плотность газа невелики. Условие на интенсивность определяется величиной пробоя в газе, моделирование которого безусловно является отдельной задачей, не связанной с теорией оптических решеток. Высокая плотность газа редко встречается в экспериментальных и численных исследованиях ввиду того, что практические задачи связаны в первую очередь с сильно разреженными газами, давление которых не превышает нескольких атмосфер. Большая область взамодействия излучения с газом, в особенности в методе ПСМ, приводит к линейному росту вычислительной сложности - например, для расчета области длиной $200\lambda$ и атмосферном давлении, требуется порядка сотни вычислительных CPU. Поэтому в задачах, в которых не учитывается обратное влияние, обычно устанавливаются периодические граничные условия на краях области, перпендикулярных направлению распространения излучения, тем самым разрешается бесконечная длина взаимодействия решетки и газа.

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

\subsection{Модель тонких оптических слоев}

\hspace{0.5cm} Решение задачи о прохождении одиночной электромагнитной волны через среду с неоднородной плотностью суть показателем преломления, даже при известном выражении $n(x)$ не имеет явного аналитического решения. Описанный ниже метод, основывающийся на дискретном представлении функции плотности газа, позволяет получить точные численное значения для прошедшей и отраженной от внутренней структуры газа волн. Приведено описание алгоритма, адаптация на случай слабоизменяющегося показателя преломления, анализ сходимости и рекомендации к применению.

\subsubsection{Распространение излучения через слоистую среду}

\hspace{0.5cm} Решение задачи об оптическом захвате газа обычно решается численно - лишь некоторые случаи (например, бесстолкновительный газ или модельный вид потенциала) имеет аналитическое решение. В настоящей работе метод ПСМ предполагает дискретное представление газовой среды в виде набора ячеек, внутри которых содержатся моделируемые пробные частиц. Одна пробная частица заменяет собой большое число реальных частиц газа. Выборка макропараметров течения производится осреднением внутри одной ячейки, и таким образом плотность газа является дискретной функцией, постоянной внутри одного вычислительного подобъема и меняющейся скачком на границе ячеек.

Такой подход, при котором функция плотности газа дискретно задается в пространстве, позволяет моделировать распространение электромагнитной волны внутри газовой среды как задачу о прохождения излучения в наборе плоских параллельных пластин с заданной диэлектрической проницаемостью. Решение этой задачи подробно приведено в \cite{Born}, и основывается на вычислении характеристических матриц сред $M_i(z)$. Для конечного набора слоев можно вычислить полные коэффициенты прохождения и отражения перемножением характеристических матриц $\prod\limits_i M_i(z_i-z_{i-1})$, а также определить локальное значение величины поля внутри каждой пластины. Данный способ, несмотря на свою универсальность и общность желательно адаптировать к настоящей задаче.

Во-первых, нас интересует прохождение волны через набор пластин, расположенных перпендикулярно направлению распространения излучения. Другими словами, волновой вектор $\vec{k}$ не имеет проекции на плоскость границы раздела пластин. В этом случае отсутствует различие между $TE$ и $TM$ типами волн. Кроме того, излучение линейно поляризовано. Во-вторых,  мы считаем, что прошедшая и отраженная волны имеют одинаковую частоту. Это не противоречит физике явления, когда два источника излучения с разными частотами формируют оптическую решетку с ненулевой скоростью, и отраженный от периодических возмущений в газе сигнал усиливает волну, распространяющуюся в противоположном направлении. Можно легко показать, что перейдя в систему отсчета, движущуюся со скоростью решетки, частоты излучения будут равны в силу эффекта Доплера. Соответственно, установить соответствующие значения частоты излучения можно уже после расчета локального распределения интенсивности в среде.

Рассмотрим одиночную границу раздела сред (рис. \ref{fig:SimpleBound}), показатели преломления которых равны соответственно $n_1$ и $n_2$. Предположим, что мы знаем все о поле в одной из сред, например, первой. Очевидно, что значение амплитуд электрической и магнитной компонент меняется (скачком) лишь на границе раздела, в то время как при прохождении однородной среды она остается неизменной.

\begin{figure}[h]
\centering
\includegraphics[width=300px]{SimpleBound}
\caption{Волны на границе раздела двух сред с разными показателями преломления.}
\label{fig:SimpleBound}
\end{figure}

Решение задачи о нахождении амплитуд поля 2 в правом полупространстве по известным амплитудам в первом (левом) полупространстве строится на условии непрерывности тангенциальных компонент полей \textbf{E} и \textbf{H} на границе $z = a$. Другими словами, ввиду отсутствия на границе раздела свободных зарядов и токов, скачок тангенциальных компонент $\{E_{\tau},H_{\tau}\}_{z = a}$ должен быть равен нулю.

Запишем условия непрерывности \eqref{fig:SimpleBound}, из которых элементарно получаются значения $E_2$ и $E'_2$.

\begin{equation}
\label{eq:SingleBound}
\begin{split}
&  E_1 e^{ik_1a} + E'_1 e^{-ik_1 a}= E_2 e^{ik_2 a}+ E'_2 e^{-ik_2 a}\\
&  n_1 (E_1 e^{ik_1 a} - E'_1 e^{-ik_1 a}) = n_2 (E_2 e^{ik_2 a} - E'_2 e^{-ik_2 a})
\end{split}
\end{equation}

Эти уравнения симметричны относительно замены индексов $1 \leftrightarrows 2$. В нашей постановке задачи мы считаем известными амплитуды с индексом 1, однако система разрешается для любых двух неизвестных. Амплитуды $E_1, E'_1$ полностью определяеют поле в всем пространстве.

Поясним, насколько данный вывод полезен для задачи о прохождения света через слоистую среду. Без ограничения общности положим, что один источник излучения находится слева от набора слоев с разными диэлектрическими проницаемостями (см. рис. \ref{fig:LayerModelE}). Если $m$ - номер последней, самой правой границы раздела, то в правом полупространстве нет волны, распространяющейся влево. Правее $m$-ой границы ничто не отражает излучение, а значит $E'_{m}(x > z_m) = 0$. Это действительно так, ведь если $E'_m$ не равна нулю, то можно подобрать параметры системы, при которых она симметрична относительно замены $z \rightarrow -z$, следовательно источник находится и справа, а мы полагали, что это не так.

\begin{figure}[h]
\centering
\includegraphics[width=450px]{LayerModelE}
\caption{Прохождение волн через слоистую среду с источником излучения, находящемся слева.}
\label{fig:LayerModelE}
\end{figure}

Пронумеруем ячейки, использующиеся в методе ПСМ, справа налево - 0-ая и $(m+1)$-ая ячейки не моделируются в процедурах ПСМ и лишь определяют граничные условия для модели оптических слоев. С физической точки зрения они задают среду, ограничивающую моделируемый газ. Ячейки, обозначенные цветом, с номерами $1,2,.. m$, содержат в себе газ постоянной плотности $N[j = 1 .. m]$ и показателем преломления $n[j] = 1 + \alpha N[j]/ \varepsilon_0$. Величины $E_{+/-}[j]$ обозначают комплексные амплитуды электрического поля с положительным и отрицательным волновыми векторами соответственно, т.е. распространяющиеся влево и вправо.

Поле $E_+[0]$ создается источником излучения и является второй известной величиной после нулевой отраженной амплитуды в правом полупространстве - амплитуда определяется интенсивностью излучения, а фазу, как будет показано ниже, мы вольны выбирать самостоятельно. Таким образом, мы обладаем достаточной информацией для определения полей во всем одномерном пространстве, и предлагается итеративно рассчитать изменение поле при переходе на каждой границе и получить комплексные величины волн, распространяющихся в положительном и отрицательном направлении во всех ячейках - слоях.

Единственной нетривиальной вещью является порядок границ для применения соотношения \eqref{eq:SingleBound} - стартовать с 0 - границы, разделяющей внешнюю среду и газ в численном подходе невозможно - неизвестна величина $E_{-}[0]$, формально определяющая коэффициент отражения от всей системы возмущенной газовой среды. Поэтому, начиная с $m$-ой границы, справа от которой известно значение $E_{-}[m+1] \equiv 0$ и определенная с точностью до константы $E_{+}[m+1]$, возможно рассчитать сначала $E_{+/-}[m]$, а затем аналогичным образом - все остальные вплоть до левого полупространства.

\textit{Замечание.} Говоря о том, что $E_{+}[m+1]$ определен с точностью до константы, мы подразумеваем следующее. Предположим, что мы знаем значения обоих полей в каждой точке пространства, и они, очевидно, определяются амплитудой и фазой волны от источника. Теперь, отнормируем каждое комплексное значение амплитуд на $E_{+}[m+1]$ - и вылетающая из среды волна (в комплексном виде) по амплитуде будет равна $(1,0)$. Чтобы вернуться к первоначальному, истинному распределению, необходимо лишь сделать нормировку на другое число, новое $\widetilde{E}_{+}[0] = {E_{+}[0] / E_{+}[m+1]}$. Поэтому в приведенном алгоритме можно установить для первичного расчета амплитуду вышедшей волны как комплексное число с нулевой фазой и единицей по модулю. Затем следует провести нормировку полученного распределения амплитуд волн на заданную граничным условием амплитуду падающего на системы излучения.

Как было сказано ранее, величины показателя преломления и апмлитуды полей внутри каждой ячейки постоянны и меняются скачком на границе. Полное выражения для поля содержит в себе также гармоническую часть $e^{i(k_jx-\omega t)}$, где $k_j = k(x_j)$ - волновой вектор волны с частотой $\omega$ в среде с показателем преломления $n(x_j)$, задающийся дисперсионным соотношением $\omega = n(x) \cdot ck$. Таким образом, полное поле $E^l_{full}$, создаваемое в исследуемой среде одним источником излучения, находящимся слева, можно записать в виде \eqref{eq:EfullSingle}.

\begin{equation}
\label{eq:EfullSingle}
E^l_{full} = E^l_{+} e^{i(k^l_jx-\omega^l t)} + E^l_{-} e^{i(-k^r_jx-\omega^r t)}
\end{equation}

Здесь мы используем обозначения $\omega^r, \omega^l$, сразу подразумевая, что падающий и отраженный лучи могут иметь разные частоты  - отражение от движущегося периодического возмущения оптической плотности приводит к изменению частоты отраженного излучения. Ввиду того, что задача об оптическом захвате самосогласованна, и скорость возмущений в каждой точке одинакова, то левый луч при отражении будет иметь частоту, равную частоте излучения от правого источника, который также формирует оптическую решетку. В случае, если в задаче отсутствует второй источник излучения, однако известна фазовая скорость возмущений, частота отраженного излучения определяется доплеровским смещением. Более подробно о прохождении встречных пучков рассмотрено ниже.

В качестве примера приведем расчет прохождения излучения через среду, показатель преломления которой меняется как гармоническая функция с заданной амплитудой (рис. \ref{main:c}). Полная длина взаимодействия $L = 200 \lambda$, одна длина волны $\lambda$ моделируется 100 ячейками, за счет чего показатель преломления изменяется очень плавно. Интенсивность прошедшей (\ref{main:a}) и отраженной (\ref{main:b}) волны уменьшаются в положительном направлении распространения излучения. Амплитуда прошедшей волны фиксирована на левой границе (единичная), в отраженной - на правой (нулевая).

\begin{figure}

\begin{minipage}{.5\linewidth}
\centering
\subfloat[]{\label{main:a}\includegraphics[width=225px]{SampleRefrEp}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\subfloat[]{\label{main:b}\includegraphics[width=225px]{SampleRefrEm}}
\end{minipage}\par\medskip
\centering
\subfloat[]{\label{main:c}\includegraphics[width=300px]{SampleRefr}}

\caption{Пространственное распределение интенсивности прошедшей и отраженной волн в среде с периодическим гармоническим возмущением показателя преломления $\Delta n_a = 3.7 \cdot 10^{-5}$.}
\label{fig:SampleRefr}
\end{figure}

\subsubsection{Пара волн в неоднородной среде}

\hspace{0.5cm} Основным применением данного метода является моделирование взаимодействия нескольких волн в сложной среде, в нашем случае - двух, формирующих поле оптической решетки. Каждый источник создает пару волн, распространяющихся в противоположные стороны - назовем их опорной (с заданной интенсивностью на входной границе) и отраженной (с заданной нулевой интенсивностью на противоположной от источника границе). Суммарное поле в среде есть суперпозиция всех четырех полей, полученных при расчете методом тонких оптических слоев. Удобно однако выделить волну, распространяющуюся в положительном направлении с амплитудой:

\begin{equation}
E_{+}(x) = E^l_{+}(x) + E^r_{+}(x) = E^l_{tr}(x) + E^r_{refr}(x)
\end{equation}

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

\begin{equation}
E_{-}(z) = E^l_{-}(z) + E^r_{-}(z) = E^l_{refr}(z) + E^r_{tr}(z)
\end{equation}

Индексы ''$tr$'' и ''$refr$'' указывают на то, что это прошедший(опорный) и отраженный сигналы, $l$ и $r$ указывают на источник излучения. Это комплексные величины со следующими граничными условиями:

\begin{equation}
\label{eq:AmplBoundCond}
\begin{split}
& E^l_{+}(x = 0) = \sqrt{I_l/2\varepsilon_0 c}\\
& E^r_{-}(x = L) = \sqrt{I_r/2\varepsilon_0 c}\\
& E^l_{-}(x = L) = E^r_{+}(z = 0) = 0
\end{split}
\end{equation}

Несмотря на то, что такая постановка задачи позволяет определить и амплитуду, и фазу поля в каждой точке пространства, эта информация избыточна. Потенциал взаимодействия частиц газа с интерференцией двух лазерных лучей с частотами $\omega_{1,2}$ - вещественная величина, пропорциональная квадрату электрического поля. В использованных обозначениях, полное электрическое поле $\mathbb{E}$ записывается как:

\begin{eqnarray*}
\mathbb{E}  & = & \frac{1}{2} \Bigl(  E_{+}(x) e^{i(k_1 x - \omega_1 t)} + E^*_{+}(x) e^{-i(k_1 x - \omega_1 t)} +\\
& &  + \  E_{-}(x) e^{i(-k_2 x - \omega_2 t)} + E^*_{-}(x) e^{-i(-k_2 x - \omega_2 t)} \Bigr)
\end{eqnarray*}

Квадрат поля $\mathbb{E}^2$, полученный напрямую из этого выражения содержит быстро осциллирующие слагаемые, пропорциональные $e^{i(\omega_1+\omega_1)t}$. Средний за некоторое разумное время (например, вычислительный шаг по времени) отклик со стороны частиц газа на такие осцилляции занулится, поэтому в дальнейшем мы будем использовать модифицированное выражение для $<\mathbb{E}^2>_\tau$, в котором высокочастотное наполнение отсутствует. При подсчете потенциала использовалось следующее выражение:

\begin{equation}
\label{eq:E2aver}
< \mathbb{E}^2 >_{\tau} = \frac{1}{2} |E_{+}(x) e^{i(k_1 x - \omega_1 t)} +  E_{-}(x) e^{i(-k_2 x - \omega_2 t)}|^2
\end{equation}

Теперь, имея удобные способы работать с получаемой методом оптических слоев информацией, необходимо верифицировать результаты на известных данных. В работе \cite{ShnKinetic} представлено аналитическое выражение для интенсивности излучения, распространяющейся в среде периодическим возмущением показателя преломления $n(x) = n_0 + \Delta n_a \cdot cos(qx-\Omega t)$. Показано, что интенсивность излучения затухает в центральной области, а также продемонстрировано превосходное согласие с численными результатами на основе модельного кинетического уравнения. На рис. \ref{fig:CompShnLayerI} изображено сравнение вида зависимости, полученной с помощью модели оптических слоев и данных \cite{ShnKinetic}.

\begin{figure}[h]
\centering
\includegraphics[width=350px]{CompShnLayerI}
\caption{Пространственное распределение интенсивности волны с положительным направлением волнового вектора в гармонически модулированной среде c амплитудой показателя преломления $\Delta n_a$ для пары лазеров одинаковой интенсивности.}
\label{fig:CompShnLayerI}
\end{figure}

Для неподвижной оптической решетки разные методы дают идентичный результат; для ненулевой скорости решетки $\xi$ зависимости отличаются незначительно - распределение интенсивности, полученное методом оптических слоев несколько ниже в центральной области. Таким образом, на основе этого сравнения можно сделать вывод, что предложенный способ расчета прохождения электромагнитных волн через оптически неоднородную среду может быть использован в задаче взаимодействия оптической решетки и газа.

\textit{Замечание.} В распределении интенсивности, полученном исследуемым методом, заметно высокочастотное наполнение (синяя и красная широкии линии на рис. \ref{fig:CompShnLayerI}). Их наличие не связано с неточностью метода или вычислительной погрешностью - оно возникает из-за самого вида показателя преломления, являющегося быстро осциллирующей на масштабе области взаимодействия функцией. Так как моель оптических слоев построена на системе уравнений Максвелла, то должен выполняться закон сохранения энергии - в частности, сохранение вектора Пойнтинга при переходе через каждую границу ячеек. Именно этот инвариант, зависящий от локальной диэлектрической проницаемости среды, дает такой эффект.

\subsubsection{Прохождение излучения через периодическую среду}

\hspace{0.5cm} Прохождение электромагнитных волн через среду с модулированным показателем преломления является основным физическим явлением, описывающим обратное воздействие среды на излучение. Отражение сигнала с длиной волны, равной длинам волн излучения, формирующих решетку, от согласованного возмущения плотности приводит к существенному затуханию проходящей волны вглубь среды. Можно показать, что брэгговское отражение от периодических структур возникает не только для гармонического вида распределения плотности, но и любой другой периодической структуры с согласованным пространственным периодом.

Рассмотрим распространение излучения от источника, находящегося слева от среды. Пара волн - опорная $E^+$ и отраженная $E^-$ будут иметь некоторое распределение амплитуд в пространстве, зависящее от длины области взаимодействия $L$, амплитуды возмущений показателя преломления $\Delta n_a$, формы возмущений и наличия случайного шума $<\delta n_a>_\lambda = 0$. Форма показателя преломления задавалась следующим образом:

1) Косинусообразный вид без учета изменения длины волны излучения на возмущениях плотности $n(x) = n_0 + \Delta n_a cos(2k_0 n_0 x)$;

2) Косинусообразный вид c учетом изменения длины волны излучения на возмущениях плотности $n(x) = n_0 + \Delta n_a cos(2k_0 n(x) x)$;

3) Меандр $n(x) = n_0 + \Delta n_a sgn(cos(2k_0 n_0 x))$;

4) Косинусообразный вид со случайным шумом $n(x) = n_0 + \Delta n_a (cos(2k_0 n_0 x) + \delta n_a(x))$.

\begin{figure}[h]
\centering
\includegraphics[width=250px]{Form_nx}
\caption{Различные виды показателя преломления с периодической структурой возмущений.}
\label{fig:Form_nx}
\end{figure}

Первые два вида используются в аналитическом решении упрощенной модели прохождения излучения через среду и практически неотличимы для небольших возмущений показателя преломления; форма со случайным шумом является аналогом пространственного распределения плотности, получаемой в расчетах методом ПСМ. Сравнение различных видов $n(x)$ для $\Delta n_a = 3.7 \cdot 10^{-5}$ и $n_0 -1 = 2.65 \cdot 10^{-4}$ представлены на рис. \ref{fig:Form_nx}.

На рис. \ref{fig:example_passE} и \ref{fig:example_refrE} представлено сравнение распределений амплитуд прошедшей и отраженной волны в средах с разной модуляцией плотности - косинусообразные возмущений без и с наложенным шумом (1 и 4) и в прямоугольной форме возмущений (3). Отметим, во-первых, что гармонический и прямоугольный вид возмущений слабо изменяет полный коэффициент отражения и прохождения волны через среду, сохраняя при этом вид зависимости амплитуд от координаты. Во-вторых, очень важно, что наличие случайного шума (синяя линия 4 на рисунке) практически не изменяет свойства прохождения излучения через среду - это позволяет заключить, что численное моделирование течения газа в оптической решетке методом ПСМ, при котором возможны случайные флуктуации плотности в расчетных ячейках и значительное локальное отклонение формы возмущений плотности от гармонических, даст корректные результаты при расчетах модели тонких оптических слоев.


\begin{figure}[!tbp]
  \centering
  \subfloat[]{\includegraphics[width=0.5\textwidth]{example_passE}\label{fig:example_passE}}
  \hfill
  \subfloat[]{\includegraphics[width=0.5\textwidth]{example_refrE}\label{fig:example_refrE}}
  \caption{Hаспределение амплитуды прошедщей (а) и отраженной (b) волн для различных видов возмущений показателя преломления.}
\end{figure}

\subsection{Численное моделирование газа в оптической решетке}

\hspace{0.5cm} Самосогласованное взаимодействие различных газов (аргон, метан) с полем оптической решетки численно моделировалось с помощью программного комплекса $Smile++$. Был проведён ряд численных экспериментов по изучению эволюции макропараметров газа (возмущений плотности, скорости, температуры) и с помощью модели тонких оптических слоев изучена эволюция распределения интенсивности излучения в газе при адиабатическом и мгновенном включении поля.

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


\subsubsection{Моделирование оптического захвата в SMILE++}

\hspace{0.5cm} Движение частиц газа при наличии внешней градиентной (дипольной) силы со стороны оптической решетки приводит к изменению их скорости. Ускорение частицы, вызванное наличием внешнего неоднородного потенциала, можно вычислить, зная градиент поля в каждой точке пространства в определенный момент времени. Предлагается использовать следующий алгоритм программной реализации:

\begin{enumerate}
\item Расчет прохождения излучения через газ с заданным распределением плотности
\item Расчет поля оптической решетки и распределения дипольной силы
\item Расчет ускорения, испытываемого частицами, в каждой точке расчетной области и изменение скорости частиц в соответствие с рассчитанным ускорением за временной шаг $\Delta t$
\item Остальные этапы метода ПСМ
\item Выборка макропараметров газа для дальнейшего расчета поля оптической решетки
\end{enumerate}

Все этапы алгоритма, за исключением п. 5, реализованы в отдельном модуле и используют общее пространство переменных комплекса $Smile++$ - данные о расчетной области, свойства газов, методы выборки макропараметров и др. 

Было создано 3 дополнительных класса: \textit{class Laser}, в котором определяются и обрабатываются характеристики двух источников излучения, \textit{class LaserField} - основной рабочий класс, в котором реализован расчет потенциала оптической решетки, и \textit{class MacroParameters} - вспомогательный класс, предназначенный для обработки макропараметров газа, необходимых для вычисления характеристик излучения в среде.

\begin{figure}[h]
\centering
\includegraphics[width=450px]{codeDiagr}
\caption{Схема программного модуля, рассчитывающего характеристики потенциала оптической решетки и его воздействия на частицы.}
\label{fig:codeDiagr}
\end{figure}

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

\textbf{1. Расчет прохождения излучения через газ с заданным распределением плотности.}

Метод расчета распространения лазерного излучения через среду с неоднородной диэлектрической проницаемостью подробно рассмотрено в первой части главы ''Результаты работы''. Именно на этом этапе последовательно определяются апмлитуды двух пар прошедших и отраженных волн с использованием информации о показателе преломления газа (обработка свойств газа осуществляется методом \textit{MacroParameters::DenstGather()} и затем \textit{MacroParameters::refractIndexGather()}).

В конце этого этапа мы получаем полную информацию о 4 волнах: распространяющихся в разные стороны от двух источников, нормированных на амплитуду входного сигнала. 

\textbf{2. Расчет поля оптической решетки и распределения дипольной силы.}

Далее, следует, во-первых, привести нормированные значения полей к реальным физическим размерностям, связанным с интенсивностью излучения (с этими параметрами работает атрибут \textit{Laser->Intensity}). Затем, используя выражение \eqref{eq:E2aver} с помощью метода \textit{LaserField::CalcInterference()} вычисляется величина $\mathbb{E}^2$, определяющая потенциал взаимодействия частиц с газом.

Сила, действующая на частицы, предполагается постоянной внутри каждой вычислительной ячейки, также как и величины поля. Это предположение оправдано всилу достаточно высокого уровня дискретизации пространственного масштаба, связанного с длиной периода оптической решетки. Сила определяется как градиент потенциала интерференционной решетки и вычисляется со вторым порядком точности $O(\Delta x^2)$ в каждой ячейке области, в том числе и на границе.

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

\textbf{3. Расчет ускорения, испытываемого частицами, в каждой точке расчетной области. Изменение скорости частиц.}

Величина изменения скорости вычислялась исходя из уравнения динамики $\ddot x = F(x,t) / m$, где сила $F(x)$ представлена в виде одномерного массива размером $NCell$ по количеству ячеек в области. Зная начальное положение частицы и ее скорость $(x,v_x | t = \tau)$, уравнение движения можно численно проинтегрировать с помощью метода Рунге-Кутта четвертого порядка за один временной шаг $\Delta t$. Шаг по времени выбирается достаточно маленьким для того, чтобы быстрые частицы не проходили сквозь потенциал оптической решетки, не замечая её периодической структуры, и испытывали плавное изменение силы вдоль траектории движения. Например, для быстрых частиц со скоростью порядка нескольких тепловых, при временном шаге $\Delta t = 10^{-12}$ сек, изменение положения вдоль $x$-координаты будет равняться $\sim 10^{-9}$ м, в то время как размер ячейки при длине волны излучения $\lambda = 532$ нм более, чем в 5 раз больше.

Несмотря на то, что количество моделируемых частиц огромно и разные частицы, находящиеся в одной и той же ячейке, могут преобретать разное изменение скорости $\Delta v_x$, основанное на их начальной скорости $v_x$, время, затраченное на главный этап работы модуля невелико, так как все значения действующей силы считаются лишь один раз и считаются постоянными на временном шаге.

\textbf{4. Остальные этапы метода ПСМ.}

После того, как вычисленное изменение скорости частицы под действием оптической решетки передается в модуль модуль \textbf{Move}, проводятся остальные, стандартные для метода ПСМ этапы моделирования. Они не зависят от наличия потенциала, поэтому полностью совпадают с уже существующими методами $Smile++$.

\textbf{5. Выборка макропараметров газа.}

После того, как все вычисления на данном временном шаге завершены, необходимо построить новое распределение показателя преломления в газе, основанное на мгновенной величине плотности частиц в каждой ячейке. На следующем шаге, обновленный массив $n[j]$ индекса рефракции будет снова использован для расчета распространения излучения через газ.

Для исследования параметров газа, например, величины возмущений плотности, на этом этапе производится вычисление амплитуды возмущения с помощью аппроксимации экспериментальных данных. В центральной части расчетной области выбирается несколько непересекающихся участков, длиной от 2 до 10 длин волн, и вычисляется апмлитуда пространственных колебаний плотности в предположении, что число частиц распределено по гармоническому закону $N(x) = N_0 + \Delta N_a \cdot cos(2kx + \phi)$, где $N_0$ - средняя плотность газа в участке, $\Delta N_a$ - искомая амплитуда возмущений, $\phi$ - фаза гармонического сигнала плотности. Амплитуда возмущений является очень важным показателем эволюции процесса, определяющим отражение излучения от среды.

\textit{Замечание.} Основным удобством использования мо оптических слоев является то, что в вычислениях не используются никакие параметры среды, кроме ее плотности (показателя преломления). В отличии от методов, описанных в \cite{ShnKinetic}, в которых делаются существенные предположения о гармоническом виде возмущений или, в случае с аналитическим решением на интенсивность излучения, о том, что амплитуда возмущения не зависит от координаты, нам нет необходимости рассчитывать $\Delta n_a$ в каждой точке пространства, что значительно повышает скорость расчета.

\subsubsection{Эволюция газа в оптической решетке}

\hspace{0.5cm} В работе большое внимание было уделено моделированию газа в неподвижной решетке - в этом случае, изначально покоящийся газ не преобретает направленной скорости и не нагревается (кроме, может быть, случая мгновенного включения излучения с высокой интенсивностью). Тогда основным меняющимся параметром газа является возмущение плотности, вызванное неоднородным внешним потенциалом.

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

В качестве среды был выбран метан (поляризуемость $\alpha_{CH4} = 2.88 \cdot 10^{-40}$ см$^2$В$^{-1}$, $m_{CH4} = 16$ а.е.м). Излучение включалось мгновенно с постоянной интенсивностью до конца расчета: $I_1(0) = I_2(L) = 10^{15} \div 10^{16}$ Вт/см$^2$. Длина области взаимодействия $L = 200\lambda$, длина волны излучения $\lambda = 532$ нм. Изначально газ находился при давлении 1 атм и температуре 293$^{\circ}$ K.

\begin{figure}[h]
\centering
\includegraphics[width=350px]{RefrIndexIcomp}
\caption{Изменение амплитуды показателя преломления для разных интенсивностей излучения.}
\label{fig:RefrIndexIcomp}
\end{figure}

На рис. \ref{fig:RefrIndexIcomp} представлены зависимости амплитуды показателя преломления $\Delta n_a \simeq \Delta N_a \alpha / 2\varepsilon_0$ для разных интенсивностей излучения. Увеличение интенсивности приводит к более сильному разделению плотности внутри газа, а следовательно и показателя преломления. На рисунке также приведено сравнение с результатом, полученным в \cite{ShnKinetic} для соответствующих условий (сиреневая пунктирная линия) - видно, что результаты имеют одинаковую ассимптотику и характерное время установления стационарного состояния, однако отличаются формой зависимости от времени. В данных, полученных в настоящей работе, отчетливо заметно наличие затухающих колебаний за время установления. Это объясняется наличием обратного влияния газа на излучение: интенсивность излучения, проходящего через среду, уменьшается в центральной области, что приводит к уменьшению силы, действующей на газ, и амлитуда возмущений снижается. Однако, при уменьшении амплитуды возмущений излучение вновь возвращается к прежднему значению, проходя через менее возмущенную среду. В таких противофазных колебаниях силы и плотности газа в области взаимодействия проявляется обратная связь между средой и излучением.

Очевидно, что распределение частиц в постоянном поле можно описать распределением Больцмана: $N(x) = N_0 \frac{exp(\alpha E^2 / 2k_BT)}{<exp(\alpha E^2 / 2k_BT)>}$, а значит можно оценить величину разделения плотности газа. Горизонтальные линии на рис. \ref{fig:RefrIndexIcomp} соответствуют величине амплитуды преломления, полученной из равновесного распределения Больцмана в заданном поле оптической решетки. Наблюдается отличное согласие оценки и численных результатов.

Ускорение и сопутствующий нагрев газа в движщейся оптической решетке можно ассимптотически описать следующим образом. Перейдем в систему отсчета, в которой решетка покоится, а газ налетает на нее с начальной скоростью $-\xi$. В пределе $t\rightarrow \infty$ газ замедляется и останавливается, взаимодействуюя с покоящимся потенциалом, в том числе засчет межмолекулярных столкновений захваченных и свободно движущихся частиц. При этом он должен нагреваться - сохраняется полная энергия системы $\frac{3}{2}k_B T_f = \frac{m\xi^2}{2}+\frac{3}{2}k_B T_0$, где $T_f, T_0$ - конечная и начальная температура газа (с 3 степенями свободы, если есть внутренние колебательные и вращательные степени свободы - то $\frac{3}{2}$ заменяется на соответствующее $\frac{i}{2}$). Таким образом, зная скорость решетки и начальную температуру газа, можно определить какую энергию получит газ за время взаимодействия. Заметим, что здесь мы не учитываем наличие нагрева засчет самого глубины потенциала - ни величина переданного импульса, ни энергии от поля газу не зависит от мощности лазерного излучения. Однако,даже при нулевой скорости решетки газ, помещенный такое поле, будет нагерваться.

На рис. \ref{fig:CompTempVx} преставлены результаты моделирования и аналитическая оценка \eqref{eq:dPdUTotal} зависимости нагрева и ускорения газа в движущейся решетке с $\xi = 411$ м/с и интенсивностью излучения $I = 10^{16}$ Вт/см$^2$. Ассимптотические оценки конечной скорости и конечной температуры газа подтверждают описанное выше предположение. Тонкая фиолетовая линия на графике температуры чуть выше зеленой (аналитическая оценка) засчет неучитываемого в модели \eqref{eq:dPdUTotal} ненулевой глубины потенциала.

\begin{figure}[h]
\centering
\includegraphics[width=400px]{CompTempVx}
\caption{Изменение амплитуды показателя преломления для разных интенсивностей излучения.}
\label{fig:CompTempVx}
\end{figure}

При получении аналитической оценки $\dot v_x$, $\dot T$ мы говорили, во-первых, об упрощении вида потенциала - представляли их в виде прямоугольных ям, а во-вторых, изменение импульса и энергии газа во времени было определено с точностью до константы. В приведенных результатах временной шаг в DSMC был в $\sqrt{2}$ раз меньше, чем в аналитической модели. Необходимо отметить, что выбор соответствия скорости изменения температуры и импульса газа между численным экспериментом и методом ПСМ может определяться универсальной константой, а значит простую аналитическую модель можно применять для быстрой оценки состояния газа в оптической решетке. В представленном случае, результаты отлично согласуются, что говорит о состоятельности такого подхода.

Сравним изменение амплитуды показателя преломления и температуры газа в неподвижной оптической решетке при разных параметрах излучения: интенсивностью $I_1 = 10^{16}$ и $I_2 = 10^{17}$  Вт/см$^2$ с мгновенным включением поля, и интенсивностью $I_3 = 10^{17}$ Вт/см$^2$ и плавным включением - физический аналог лазерного импульса с длительностью $5$ нс и центром импульса в $5$ нс спустя начало эксперимента. В последнем случае, зависимость интенсивности обоих лазеров от времени имеет гауссовский профиль, а следовательно и сила, действующая на частицы, будет увеличиваться до момента $5$ нс и уменьшаться после.



\subsubsection{Пространственное распределение и эволюция излучения}

Эволюция распределения излучения в возмущенной среде зависит от множества параметров - свойств газа, интенсивности самого излучения, размера области взаимодействия. Несмотря на то, что полученные в ранних работах результаты затухания интенсивности в центральную область говорят о ''тонкости'' эффектов - отличие действующей силы на газ с и без учета обратного влияния газа составляет порядка $1 \%$ на масштабах взаимодействия $100-1000 \lambda$, изучение того, как затухает излучение вглубь среды крайне важно в макроспопических задачах с $L \sim 1$ см. Кроме того, отличается ли эволюция распределения поля внутри газа при мгновенном и адиабатическом включении потенциала?

Рассмотрим изменение амплитуд пары прошедшей и отраженной волн, создаваемой одним источником, например, левым. Анализ эволюции распределения величин $E^+_l(x,t)$ и $E^-_l(x,t)$ позволит более детально описать процессы, связанные с установлением равновесия в системе. На рис. \ref{fig:Ep_comp} представлены распределения нормированной амплитуды поля $E_l^+$ - опорного луча от левого источника. В нулевой момент времени, когда среда не возмущена, волна проходит через среду без отражения и амплитуда равна единице везде. Также, из условия нормировки, амплитуда равна единице на входной границе области $x = 0$ во все время взаимодействия. 

\begin{figure}

\begin{minipage}{.5\linewidth}
\centering
\subfloat[]{\label{Ep_comp:a}\includegraphics[width=250px]{Ep_2000L_1016W_Ar}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\subfloat[]{\label{Ep_comp:b}\includegraphics[width=250px]{Ep_200L_1017W_CH4}}
\end{minipage}\par\medskip
\centering
\subfloat[]{\label{Ep_comp:c}\includegraphics[width=250px]{Ep_2000L_1016W_Ar_gaus}}

\caption{Пространственное распределение нормированной амплитуды прошедшей волны. a) $L = 2000\lambda, I = 10^{16}$ Вт/см$^2$, b) $L = 200\lambda, I = 10^{17}$ Вт/см$^2$, c) $L = 2000\lambda, I = 10^{16}$ Вт/см$^2$, длительность импульса 5 нс.}
\label{fig:Ep_comp}
\end{figure}

При низкой интенсивности $I(x=0) = 10^{16}$ Вт/см$^2$ (рис. \ref{Ep_comp:a}) амплитуда достаточно быстро приходит к почти стационарному распределению, даже при большом размере области $L = 2000\lambda$. Характерное время колебаний амлитуды $\tau_E$ составляет 0.2 нс, что соответствует характерному времени между столкновениями частиц при данных условиях $\tau_{col} = l_{col}/v_T \simeq 0.15$ нс. На рис.\ref{Ep_comp:b} представлен другой случай с $I(x=0) = 10^{17}$ Вт/см$^2$ и $L = 200\lambda$. Более высокая интенсивность излучения приводит к значительно большей амплитуде возмущений плотности, и амплитуда осциллирует с меньшим дикрементом затухания. В отличие от случая $a$ с меньшей интенсивностью, при которой можно определить наличие 1-2 периодов колебания, в случае $b$ периодов колебаний не менее 4. Тем не менее, время установления стационарного распределения амплитуды одинакова во обоих случаях и составляет $\tau_{rel} = 1 - 1.5$ нс. Это время совпадает с другим временным масштабом задачи - временем колебаний молекул между стенками оптической решетки с тепловой скоростью $\tau_{\lambda} \sim \frac{\lambda}{v_T} = 1.3$ нс.


\begin{figure}

\begin{minipage}{.5\linewidth}
\centering
\subfloat[]{\label{main:a}\includegraphics[width=240px]{Em_2000L_1016W_Ar}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\subfloat[]{\label{main:b}\includegraphics[width=240px]{Em_200L_1017W_CH4}}
\end{minipage}\par\medskip
\centering
\subfloat[]{\label{main:c}\includegraphics[width=250px]{Em_2000L_1016W_Ar_gaus}}

\caption{Пространственное распределение нормированной амплитуды отраженной волны. a) $L = 2000\lambda, I = 10^{16}$ Вт/см$^2$, b) $L = 200\lambda, I = 10^{17}$ Вт/см$^2$, c) $L = 2000\lambda, I = 10^{16}$ Вт/см$^2$, длительность импульса 5 нс.}
\label{fig:Em_comp}
\end{figure}

В случае, если излучение включается медленно (длительность гауссовского импульса 5 нс, центр импульса 3 нс с момента $t = 0$), амплитуда показателя преломления, как было показано ранее, плавно изменяется, и при времени, большем продолжительности импульса, возмущения плотности исчезают. Соответствующим образом ведет и амплитуда прошедшей волны (рис. \ref{Ep_comp:c}) - максимальное ослабление достигается в момент времени максимального возмущения и максимальной интенсивности. Заметим, что гауссовский профиль $E^+_l(x)$ определяется свойствами среды, а не модуляцией энергии лазера из-за формы импульса, так как нормировка проводится на входную интенсивность излучения в соответствующий момент времени.

Рис. \ref{fig:Em_comp} содержит распределение амплитуды отраженной волны $E^-_l(x)$, создаваемую левым источником. Анализ вида зависимостей аналогичен предшествующему для $E^+_l(x)$.

\begin{figure}[t]
\begin{multicols}{2}
\hfill
\includegraphics[width=83mm]{I1_200L_1016W_CH4}
\hfill
\caption{Интенсивность излучения в газе метана при длине области взаимодействия $L = 200\lambda$. $I(x = 0) = 10^{16}$ Вт/см$^2$.}
\label{fig:I1_200L_1016W_CH4}
\hfill
\includegraphics[width=83mm]{I1_2000L_1016W_Ar}
\hfill
\caption{Интенсивность излучения в газе аргона при длине области взаимодействия $L = 2000\lambda$. $I(x = 0) = 10^{16}$ Вт/см$^2$.}
\label{fig:I1_2000L_1016W_Ar}
\end{multicols}
\end{figure}

Аналогично эволюции интенсивности, интерес представляет рассмотрение эволюции интенсивности волны, распространяющейся в положительном направлении $I_l(x,t)$ - она складывается из опорного (прошедшего) луча, создаваемого левым источником и отраженного луча от правого. Можно предположить, что при увеличении области взаимодействия интенсивность должна изменяться сильнее в центральной части газа, так как отражение будет происходить от большего числа возмущений. На рис.\ref{fig:I1_200L_1016W_CH4} и \ref{fig:I1_2000L_1016W_Ar} представлено пространственное распределение интенсивности одной волны в разные моменты времени для излучения мощностью ${10^{16}  W/cm^2}$ при разных размерах области взаимодействия $L_1 = 200 \lambda$, $L_2 = 2000 \lambda$. Интенсивность нормирована на входное значение $I_l(x=0)$. В обоих случаях, она быстро уменьшается за время $\tau_I \simeq 0.2$ нс на $0.15\%$ и $2\%$ соответственно, а затем возвращается к почти номинальному значению и приходит к динамическому стационару. Заметим, что распределение интенсивности излучения не выходит на полностью стационарное решение - это связано с наличием описанной выше обратной связи между распределением плотности и отражением от нее волн. Однако, при достаточной величине столкновительной диссипации в газе (когда диффузия достаточно велика), при больших временах взаимодействия возможно установление квазистационара, при котором небольшие колебания распределения интенсивности и амплитуды возмущения плотности будут малозаметны.


Обсуждаемые колебания распределения интенсивности должны быть более выражены, если увеличить энергию излучения - это приведет к большему разделению плотности газа в потенциале, а значит и усилению отражения от периодических возмущений показателя преломления. На рис. \ref{fig:I1_200L_1017W_CH4} предстален расчет для $L = 200\lambda$ и интенсивностью в 10 раз выше, чем в предыдущих случаях ($I_1(0) = I_2(L) = 10^{17}$ Вт/см$^2$).

\begin{figure}[h]
\centering
\includegraphics[width=400px]{I1_200L_1017W_CH4}
\caption{Интенсивность излучения в газе метана при длине области взаимодействия $L = 200\lambda$. $I(x = 0) = 10^{17}$ Вт/см$^2$.}
\label{fig:I1_200L_1017W_CH4}
\end{figure}

Даже при достаточно небольшом размере области, интенсивность изменяется более чем на $2\%$ и характерное время колебаний составляет $\tau_I \sim 0.3$ нс. В работе \cite{ShnKinetic} сделано предположение, что в описываемой системе можно ожидать периодический затухающий процесс колебаний плотности суть показателя преломления с характерным временем порядка среднего времени между столкновениями в газе $\tau_c$. В исследуемом случае, это время составляет примерно $0.2\div0.4$ нс, что отлично согласуется с временем колебаний интенсивности $\tau_I$.