\documentclass[twocolumn]{jetpl}
\usepackage[russian]{babel}
% \bibliographystyle{maik}
\usepackage[koi8-r]{inputenc}
\usepackage[T1,T2A]{fontenc}
\usepackage{cite}
% \usepackage[sort]{cite}
% \usepackage{multicol}
\usepackage{physics}
\usepackage{mathtools}
\usepackage{xcolor}


% \usepackage{biblatex}
%%% article in English
\rus

%%% declaration of a new mathematical operator
\DeclareMathOperator{\sign}{sign}

%%% article title
\title{Нелинейный отклик разреженных газов на ультрафиолетовый фемтосекундный импульс}

%%% article title - for colontitle (at the top of the page)
\rtitle{Нелинейный отклик разреженных газов на ультрафиолетовый фемтосекундный импульс}

%%% article title - for table of contents (usualy identical with \title)
\sodtitle{Нелинейный отклик разреженных газов на ультрафиолетовый фемтосекундный импульс}

%%% author(s) ( + e-mail)
\author{
Н.\,Р.\,Врублевская$^{+}$, Д.\,Е.\,Шипило$^{+*}$, И.\,А.\,Николаева$^{+*}$, Н.\,А.\,Панов$^{+*}$, О.\,Г.\,Косарева$^{+*}$\thanks{e-mail: kosareva@physics.msu.ru}
}

%%% author(s) - for colontitle (at the top of the page)
\rauthor{Врублевская, Шипило, Николаева, Панов, Косарева}

%%% author(s) - for table of contents
\sodauthor{Н.\,Р.\,Врублевская$^{+}$, Д.\,Е.\,Шипило$^{+*}$, И.\,А.\,Николаева$^{+*}$, Н.\,А.\,Панов$^{+*}$, О.\,Г.\,Косарева$^{+*}$\thanks{e-mail: kosareva@physics.msu.ru}}

%%% author's address(es)
\address{~\\$^+$Физический факультет МГУ имени~М.\,В.~Ломоносова, Москва, Ленинские горы, 1, стр. 62, 119991, Россия\\~\\
$^*$Физический институт им. П.\,Н.~Лебедева РАН, Москва, Ленинский пр-т, 53, 119991, Россия}

%%% dates of submition & resubmition (if submitted once, second argument is *)
\dates{2 февраля 2023 г.}{\today}

%%% abstract
\abstract {Квантовомеханические расчеты нелинейного отклика одномерной квантовой системы, воспроизводящей энергетическую структуру ксенона, на ультрафиолетовый фемтосекундный импульс с интенсивностью 1--100\,ТВт/см$^2$ показали дисперсию коэффициента кубической нелинейности в диапазоне 266--400\,нм и его зависимость от интенсивности, исключающую описание отклика связанных электронов в~виде $\chi^{(3)}E^3$. 
Вычисление поляризации на базе одномерной квантовой модели может быть использовано при моделировании распространения ультрафиолетового фемтосекундного излучения в~газе.
} 

%%% PACS numbers
\PACS{42.65.-k, 42.65.An}


\begin{document}

\maketitle

\textbf{1. Введение.} К настоящему времени проведено детальное экспериментальное и теоретическое исследование  распространения высокоинтенсивного фемтосекундного лазерного излучения инфракрасного диапазона в газах как при филаментации в объеме газовой среды~\cite{ReviewKand09} так и в газонаполненном волокне~\cite{Zheltikov2004fiber_book}.
Для теоретического описания нелинейного отклика газовых сред на фемтосекундный импульс инфракрасного диапазона применяется \textit{феноменологический} подход~\cite{couairon2011practitioner}, основанный на представлении нелинейности в виде суперпозиции тока электронов, освободившихся в актах многофотонной/туннельной ионизации~\cite{brunel1990harmonic}, и~производной по времени $t$ нелинейной поляризации среды $P_{nl}$, соответсвующей ангармоничному движению связанных электронов.
Поляризация $P_{nl}$, в основном, определяется мгновенным кубическим по полю $E$ вкладом $P_{nl} =\chi^{(3)} E^3$~\cite{kosareva2011arrest,polynkin2011experimental,wahlstrand2011optical,volkova2012polarisation}, коэффициент нелинейности третьего порядка \mbox{$n_2 \propto \chi^{(3)}$} определяется экспериментально~\cite{zahedpour2015measurement,Chin05,kompanets2020nonlinear}.
Для расчетов распространения в~молекулярных газах модель дополняют описанием резонансного линейного~\cite{PanovPRA16,Geints10mkm2017,panov2019nonlinear,kompanets2020nonlinear} и инерционного кубического~\cite{PlatonenkoLP93,NibberingInertialJOSAB97,PalastroRaman2012,panov2019nonlinear,kompanets2022} откликов.

Строго говоря, феноменологический подход к~описанию нелинейности не обоснован математически, его применимость обусловлена существенно различным числом фотонов, обеспечивающих отклик связанных электронов ($3\hbar\omega$) и ионизацию среды ($8\hbar\omega$ для молекулы О$_2$ при центральной длине волны излучения $\lambda = 2\pi c/\omega = 800$\,нм, $c$~--- скорость света).
Для фемтосекундных импульсов с центральной длиной волны в ближнем и среднем инфракрасном диапазоне феноменологический подход воспроизводит экспериментальные результаты.
Например, для наиболее часто используемого в экспериментах излучения с центральной длиной волны около 800\,нм (лазер на титан-сапфире) оцененная из измерений и моделирования пиковая интенсивность в~филаменте в~воздухе составляет около 100\,ТВт/см$^2$~\cite{Chin00,Kosareva09,xu2012intensity_by_fluorescence,Reevaluation15}. 

В плазменном канале ультрафиолетового филамента возможно усиление радиочастотного и~терагерцового излучения благодаря относительно узкому спектру фотоэлектронов~\cite{bogatskaya2020jetp}.
Однако параметры филаментов ультрафиолетового диапазона не столь детально изучены, как для инфракрасного излучения.
В частности, неизвестен даже порядок интенсивности импульса на длине волны $\sim$250\,нм (третья гармоника лазера на титан-сапфире) в~филаменте: согласно работе~\cite{TzorUV:OL00} она составляет около 0.1\,ТВт/см$^2$, тогда как согласно~\cite{PhysRevLett.88.135003}~--- 20\,ТВт/см$^2$.
Вызывает сомнения сама возможность применения феноменологического подхода~\cite{PhysRevLett.88.135003,fedorov2008wavelength_scan,ShipiloUV2017,shutov2019humidity_plasma} для теоретического описания взаимодействия высокоинтенсивного фемтосекундного лазерного излучения ультрафиолетового диапазона с газовыми средами.
Действительно, для излучения с~длиной волны $\sim$250\,нм ионизация таких газов как кислород или ксенон становится трехфотонной~\cite{fedorov2016keldysh}.
Поэтому отклик электронов как связанных в атомах или молекулах, так и освобожденных в актах ионизации~\cite{lewenstein1994theory,brunel1990harmonic} становится кубичным по полю.

Косвенное подтверждение ограниченной применимости феноменологического подхода к описанию нелинейности газов в ультрафиолетовом диапазоне представлено в работе~\cite{CouaironMIRn2}, в которой построена аппроксимация зависимости коэффициента $n_2$, полученного на основе квантовых вычислений, от длины волны $\lambda$ формулой селлмейеровского типа.
Сингулярность в аппроксимации $n_2(\lambda)$ достигалась на длине волны, соответствующей трети потенциала ионизации,~--- то есть когда нелинейная ионизация становится трехфотонной.
Тем самым, для фемтосекундного излучения ультрафиолетового диапазона с~частотой вблизи трети потенциала ионизации само определение $n_2$ может терять физический смысл.

%\begin{figure}
%    \centering
%    \includegraphics[width=\linewidth]{Fig1.eps}
%    \caption{Рис. 1. (Слева) Потенциал $U(x)$ и волновые функции связанных состояний $\Psi_j(x)$. (Справа)~Атомные уровни ксенона. 
%    }
%    \label{fig:pit}
%\end{figure}

Стоит отметить, что квантовомеханическое описание нелинейного отклика полупроводников и~диэлектриков обычно базируется на системе уравнений Максвелла---Блоха, см.,~например,~\cite{pfeiffer2020Bloch,gaarde2020structure}. 
В случае разреженного одноатомного газа, однако, отсутствует необходимость учитывать зонную структуру, анизотропию среды, взаимодействие с~резервуаром и~термализацию плазмы, 
поэтому можно отказаться от такого формализма в~пользу нестационарного уравнения Шредингера~\cite{volkova1992dynamics,stremoukhov2022quasi}.

В настоящей работе мы, используя одномерную квантовомеханическую модель~\cite{volkova1994tunneling} взаимодействия света с потенциальной ямой с уровнями энергии, приближенно соответствующими основному и~возбужденным состояниям атома ксенона, исследуем нелинейный отклик, наводимый в среде фемтосекундным импульсом.
В~широком диапазоне интенсивностей импульса от $0.1$ до $100$\,ТВт/см$^2$ и~его центральных длин волн от 266 до 1500\,нм из численного решения нестационарного уравнения Шредингера нами получен нелинейный отклик среды, состоящей из таких невзаимодействующих между собой квантовых систем. Установлено, что в~инфракрасном диапазоне найденная нами нелинейная поляризация соответствует отклику, определенному из феноменологического подхода.
Для ультрафиолетовых фемтосекундных импульсов это согласие нарушается: уже при низкой интенсивности 2\,ТВт/см$^2$ нелинейный отклик не может быть аппроксимирован кубом поля, а~с~ростом интенсивности до 25\,ТВт/см$^2$ рассчитанные квантовомеханически и феноменологически зависимости нелинейной поляризации от времени осциллируют со сдвигом фазы в~четверть оптического периода и~имеют различный частотный состав. 

\textbf{2. Одномерная квантовомеханическая модель взаимодействия фемтосекундного импульса с веществом.}
Пусть $U(x)$~--- одномерная потенциальная яма.
Для определения связанных состояний потенциала $U(x)$ мы использовали итерационный алгоритм, подробно описанный в Приложении. 
Вариацией параметров мы построили такой потенциал (здесь и далее в формулах использованы атомные единицы)
\begin{equation}
U(x) = -\frac{0.5625}{\sqrt{x^2 + 0.63^2}} \exp \Big[-\Big(\frac{x}{8}\Big)^{16}\Big],
\label{eq:potential}    
\end{equation}
что его связанные состояния $\ket{\Psi_j}$, $j=0,1,2$ с~энергиями $W_0 = -12.08$\,эВ, $W_1 = -2.93$\,эВ и $W_2 = -1.17$\,эВ воспроизводят энергетическую структуру атома ксенона (см. рис. 1). 

Для описания взаимодействия света с веществом использовалось нестационарное уравнение Шредингера для волновой функции $\Psi(x,t)$ с начальными условиями $\Psi(x,t \xrightarrow{}-\infty) = \ket{\Psi_0}$: 
\begin{equation}
i\frac{\partial\Psi}{\partial t} = 
-\frac{1}{2} \frac{\partial^2 \Psi}{\partial x^2} + U(x)\Psi - E(t)x\Psi,
\label{eq:shred}    
\end{equation}
где $ E(t) = \sqrt{I_0} \exp(-t^2/[2\tau_0^2])\sin{\omega_0 t}$~--- электрическое поле лазерного импульса длительностью $\tau_0 = 5$\,фс.
Интенсивность излучения $I_0$ варьировалась от $10^{11}$ до $10^{14}$\,Вт/см$^2$ (от 3$\times 10^{-6}$ до 3$\times 10^{-3}$ атомной интенсивности), частота $\omega_0$ соответствовала длинам волн от 1500 до 266~нм, то есть от низкочастотного туннельного предела до трехфотонной ионизации.

Расчёты проводились на видеокарте NVIDIA GeForce RTX 3080 с использованием технологии CUDA. Общий размер временной области составлял $20\,\tau_0 = 100$\,фс, пространственной области~--- 8192\,aт.\,eд.
Количество узлов временной и~пространственной сеток $N = 2^{16}$ обеспечивало разрешение по времени $\Delta t = 1.5$\,ас и по пространству $\Delta x = 0.125$\,ат.\,ед. Характерное время расчета составляло около 5\,минут. 

Полученная при численном интегрировании уравнения~\eqref{eq:shred} вероятность ионизации  $\eta(t)=1-\sum_j|\braket{\Psi_j}{\Psi}|^2$ находится в согласии с~результатами расчета по формуле Переломова---Попова---Терентьева. Для центральных длин волн 266, 800 и 1500~нм зависимости вероятности перехода из связанного состояния одномерной квантовой системы в свободное от интенсивности приведены на рис.~2. 
Они соответствуют многофотонному пределу для ультрафиолетовых импульсов и туннельному для импульсов ближнего и среднего инфракрасного диапазона. 

%\begin{figure}
%    \centering
%    \includegraphics[width=\linewidth]{Fig2.eps}
%    \caption{Рис. 2. Вероятности ионизации в широком диапазоне значений пиковой интенсивности импульса длительностью 5\,фс (символы) и её аппроксимация в многофотонном и туннельном пределах (кривые).
%    }
%    \label{fig:ionization}
%\end{figure}

\textbf{3. Нелинейный отклик одномерной квантовой системы.}
Поскольку мы рассматриваем разреженные газовые среды, 
макроскопическая поляризация $P$ равна произведению дипольного момента атома газа и~концентрации частиц. Произведем нормировку макроскопической поляризации на концентрацию частиц. Нормированная поляризация будет равна дипольному моменту нашей квантовой системы:
\begin{equation}
P(t) = -\braket{\Psi|\hat{x}}{\Psi}. 
\label{eq:polarization_q}    
\end{equation}
Поляризация $P(t)$ содержит как линейный по полю, так и нелинейный вклады связанных электронов, а~также вклад электронов континуума и интерференционные члены, разделить которые в~общем случае невозможно. 
Определим эффективный коэффициент кубической нелинейности $n_2\propto\chi^{(3)}$ из наилучшей аппроксимации поляризации~\eqref{eq:polarization_q} зависимостью 
\begin{equation}
P(t) = \chi^{(1)}(t)\otimes E(t) + \chi^{(3)}E^3(t),
\label{eq:polarization_c}
\end{equation}
где знаком $\otimes$ обозначена свертка. Используя для аппроксимации расчеты при нескольких значениях пиковой интенсивности $I_0$, можно выделить нелинейную часть поляризации. При этом дисперсия линейного отклика газовой среды учитывается для всего диапазона пиковой интенсивности.

По найденным из аппроксимации~\eqref{eq:polarization_c} значениям $\chi^{(3)}$ на рис.~3 построены зависимости эффективного коэффициента кубической нелинейности $n_2$ как функции интенсивности для трех длин волн $\lambda=1500$, 800 и 266\,нм.
Эти значения примерно на порядок меньше известных из экспериментов~\cite{zahedpour2015measurement,Chin05,kompanets2020nonlinear}  значений коэффициента керровской нелинейности газов в инфракрасном диапазоне $n_2 \sim 10^{-19}$\,см$^2$/Вт. 
Для импульсов с $\lambda=1500$, 800\,нм эффективный коэффициент $n_2$ постоянен вплоть до интенсивностей 30--40\,ТВт/см$^2$, после чего становится отрицательным из-за существенной доли ионизованных атомов ($\eta \gtrsim 10^{-3}$, см.~рис.~2).
В~случае $\lambda=266$~нм рост ионизации с интенсивностью происходит <<постепенно>>, и~уже на интенсивностях более 5~ТВт/см$^2$ нельзя считать $n_2$ постоянным. 

%\begin{figure}
%    \centering
%    \includegraphics[width=\linewidth]{Fig3.eps}
%    \caption{Рис. 3. Зависимости эффективного коэффициента кубической нелинейности от интенсивности для различных длин волн. Для инфракрасных импульсов зависимость можно считать постоянной, для ультрафиолетового~--- линейно убывающей (линии).  
%    }
%    \label{fig:n2}
%\end{figure}

%\begin{figure}
%    \centering
%    \includegraphics[width=\linewidth]{FigDisp.eps}
%    \caption{Рис. 4. Зависимость эффективного коэффициента кубической нелинейности, оцененного согласно аппроксимации~\eqref{eq:polarization_c} из квантовомеханических расчетов~\eqref{eq:shred}, от центральной частоты фемтосекундного импульса с~пиковой интенсивностью 10\,ТВт/см$^2$ (черные точки) и~ее аппроксимация формулой селлмейеровского типа~\eqref{eq:sellmeier} (красная линия).
%    }
%    \label{fig:n2_on_freq}
%\end{figure}

Для демонстрации дисперсии кубической нелинейности на рис.~4 построена зависимость эффективного коэффициента $n_2$ от центральной частоты лазерного импульса. 
В~согласии с данными моделирования~\cite{CouaironMIRn2} она может быть аппроксимирована формулой селлмейеровского типа~\cite{fedorov2008wavelength_scan}, если в~качестве резонансных частот взять 1/3 от резонансных частот линейной дисперсии, которые соответствуют переходам из основного состояния в~первое возбужденное и~в~континуум:
\begin{equation}
n_2(\omega_0) = \frac{A}{\Omega_{A}^2 - \omega_0^2} + \frac{B}{\Omega_{B}^2 - \omega_0^2},
\label{eq:sellmeier}
\end{equation}
где $\Omega_A = |W_0|/3=4.03$\,эВ, $\Omega_B = |W_0-W_1|/3=3.05$\,эВ.
На частотах выше 600--700\,ТГц (длинах волн меньше 430--500\,нм) значения $n_2$ существенно меняются с~изменением частоты в~пределах спектральной ширины фемтосекундного импульса и~становятся отрицательными. 
Аппроксимация~\eqref{eq:polarization_c} перестает быть физически осмысленной, поскольку нелинейный отклик третьего порядка является запаздывающим и/или зависит от интенсивности излучения.  

%\begin{figure}
%    \centering
%    \includegraphics[width=\linewidth]{1500nm.eps}
%    \caption{Рис. 5. Нелинейная поляризация атомной системы, полученная в~квантовомеханическом расчете (красные линии) и~на основе феноменологической модели (черные линии) при воздействии импульса на длине волны 1500\,нм с~интенсивностью (a) 15 и (b,c) 59\,ТВт/см$^2$. 
%    (d) Поле лазерного импульса.
%    }
%    \label{fig:1500nm}
%%\end{figure}

%\begin{figure}
%    \centering
%    \includegraphics[width=\linewidth]{266nm.eps}
%    \caption{Рис. 6. То же, что на рис. 5, для импульсов на длине волны 266\,нм с~интенсивностью (a)~2 и (b,c)~25\,ТВт/см$^2$. 
%    }
%    \label{fig:266nm}
%\end{figure}

Для инфракрасных импульсов [1500\,нм, рис.~5, рис.~7(a)] нелинейную поляризацию квантовой системы можно с~хорошей точностью воспроизвести в~феноменологической модели, подобрав коэффициент кубической нелинейности и~предэкспоненциальный фактор в скорости ионизации.
Высокочастотные осцилляции поляризации на заднем фронте импульса соответствуют рекомбинационной генерации гармоник~\cite{serebryannikov2014quantum,zheltikov2021ufn_ru}.
Для ультрафиолетового импульса (266\,нм, рис.~6) 
можно добиться совпадения амплитуд $P(t)$, рассчитанных двумя способами, но при этом у них существенно отличается фаза, внутрипериодная динамика и~частотный спектр [рис.~7(b)].
Даже при низкой интенсивности 2\,ТВт/см$^2$ [рис.~6(a)], когда доля ионизованных атомов $\eta\sim 10^{-5}$, в~квантовой системе 
поляризация не пропорциональна кубу электрического поля и запаздывает относительно него.
С~ростом интенсивности до 25\,ТВт/см$^2$ запаздывание нелинейного отклика, найденного при численном решении уравнения Шредингера, относительно отклика, расчитанного в~рамках феноменологического подхода, увеличивается, достигая четверти оптического периода [рис.~6(b),~6(c)].

%\begin{figure}

%    \centering
 %   \includegraphics[width=\linewidth]{Fig7.eps}
  %  \caption{Рис. 7. Спектры нелинейной поляризации для (a) инфракрасного и (b) ультрафиолетового импульсов. Спектры получены из зависимостей, построенных на рис. 5(c) и 6(c), соответственно.
   % Серые вертикальные линии указывают центральные частоты импульсов.
    %}
 %   \label{fig:spectra}
%\end{figure}

\textbf{4. Заключение.} На основе разработанной одномерной квантовомеханической модели взаимодействия лазерного излучения с~веществом исследован нелинейный отклик атома на фемтосекундное излучение с~центральными длинами волн от 266 до 1500\,нм.
В ультрафиолетовой части спектра эффективный коэффициент кубической нелинейности демонстрирует дисперсию и~сильную зависимость от интенсивности, то есть $\chi^{(3)}$ не может считаться константой.
Даже при интенсивности $\sim$1\,ТВт/см$^2$ нелинейный отклик на ультрафиолетовый импульс не является кубическим по полю.
С ростом интенсивности до $\sim$20\,ТВт/см$^2$ отклонение результатов квантовомеханического расчета нелинейной поляризации от результатов расчета в рамках феноменологического подхода увеличивается.

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

\textbf{Благодарности.}
Работа поддержана грантом Российского научного фонда (21-49-00023). 
Работа Д.\,Е.\,Шипило поддержана стипендией Президента РФ молодым ученым и аспирантам (СП-3450.2022.2).
Работа И.\,А.\,Николаевой и Н.\,Р.\,Врублевской поддержана стипендиями Фонда развития теоретической физики и математики <<БАЗИС>> \mbox{(21-2-10-55-1 и 22-2-1-41-1)}.

\textbf{Приложение: итерационный алгоритм определения связанных состояний потенциальной ямы.}
Для опеределения уровней энергии заданного потенциала $U(x)$ и волновых функций $\ket{\Psi_j}$
к~произвольно выбранному начальному приближению $\ket{\Psi_0}^{(0)}$ многократно применялся оператор \mbox{$(1 - \alpha\hat{H})$}:
\begin{equation}
\ket{\Psi_0}^{(k + 1)} = \hat{N}(1 - \alpha\hat{H})\ket{\Psi_0}^{(k)},    
\end{equation}
где $\hat{H} = -\frac{1}{2}\frac{\partial^2}{\partial x^2} + U$~--- гамильтониан, $\hat{N}\ket{\varphi}=\ket{\varphi}/\sqrt{\braket{\varphi}{\varphi}}$~--- нормировка волновой функции, а параметр $\alpha \ll 1$ подбирался таким образом, чтобы самые высокие состояния континуума на заданной сетке затухали. 
Поскольку самое высокое состояние континуума соответствует частоте Найквиста, достаточно  потребовать $\alpha < (2\Delta x/\pi)^2$, где $\Delta x$~--- шаг сетки.
Стоит отметить, что требовалось достаточно большое количество итераций, до 100\,000, чтобы полученная таким образом функция основного состояния $\ket{\Psi_0}$ была ортогональна остальным волновым функциям и не приводила к артефактам в численном решении нестационарного уравнения~\eqref{eq:shred}. Высшие связанные состояния  $\ket{\Psi_1}$ и $\ket{\Psi_2}$ находились аналогично, с той разницей, что в итерационный алгоритм было добавлено условие ортогональности найденным ранее состояниям:
\begin{equation}
\ket{\Psi_1}^{(k + 1)} =  \hat{N}
\Big(1 - \ket{\Psi_0}\bra{\Psi_0}\Big)
\left(1 - \alpha\hat{H}\right) \ket{\Psi_1}^{(k)}
\end{equation}
\begin{equation}
\ket{\Psi_2}^{(k + 1)} =  \hat{N}
\left(1 - \sum_{j=0}^1\ket{\Psi_j}\bra{\Psi_j}\right)
\left(1 - \alpha\hat{H}\right) \ket{\Psi_2}^{(k)}    
\end{equation}


\bibliographystyle{ieeetr}
\bibliography{filamentation.bib}

\end{document}