[a,b]Vsevolod Chestnov

#  Recent progress in intersection theory for Feynman integrals decomposition

###### Abstract

High precision calculations in perturbative QFT often require evaluation of
big collection of Feynman integrals. Complexity of this task can be greatly
reduced via the usage of linear identities among Feynman integrals. Based on
mathematical theory of intersection numbers, recently a new method for
derivation of such identities and decomposition of Feynman integrals was
introduced and applied to many non-trivial examples.

In this note based on [1] we discuss the latest developments in algorithms for
the evaluation of intersection numbers, and their application to the reduction
of Feynman integrals.

##  1 Introduction

Feynman integrals (FIs) are a corner stone of the perturbative Quantum Field
Theory (pQFT), at least in its contemporary formulation. Within that branch of
research one fruitful discovery were the integration-by-parts identities
(IBPs) [2, 3], i.e. linear identities among FIs. For a given pQFT problem
(such as, for example, scattering amplitudes) IBPs allow to reduce an infinite
set of contributing FIs to a linear combination of finite number of basic
objects known as the master integrals (MIs).

Recently a novel framework [4, 5, 6, 7, 8] based on the twisted cohomology
theory [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19] was proposed to describe
relations among FIs. It was shown that (for fixed topology of the
corresponding Feynman graphs) FIs form a finite dimensional vector space
endowed with a scalar product called the intersection number. Among other
applications, this structure then helped to derive novel algorithms for the
direct projection of FIs onto the basis of MIs.

In Section 2 review the basics of twisted cohomology theory and computation of
intersection numbers. Then in Section 3 we present another algorithm111 This
section is based on the joint work with Federico Gasparotto, Manoj K. Mandal,
Pierpaolo Mastrolia, Saiei J. Matsubara-Heo, Henrik J. Munch, Nobuki Takayama.
[1] for reduction of FIs exploiting the connection with the Gel’fand-Kapranov-
Zelevinsky (GKZ) hypergeometric systems and the secondary equation [19].

##  2 Twisted cohomology

Here we review some aspects of the twisted cohomology and intersection theory
see also [20, 21, 22, 23] and [24]. Our central subject of study is going to
be generalized hypergeometric integrals of the form:

$\displaystyle\mathcal{I}=\int_{\mathcal{C}}u(x)\,\varphi(x)\ ,$ (1)

where $u(x)=\prod_{i}\mathcal{B}_{i}^{\gamma_{i}}$ is a multivalued function,
$\mathcal{C}$ is an $n$-dimensional integration contour such that
$\prod_{i}\mathcal{B}_{i}(\partial\mathcal{C})=0$, and
$\varphi\equiv\hat{\varphi}\>\mathrm{d}x_{1}\wedge\ldots\wedge\mathrm{d}x_{n}$
is a holomorphic $n$-form (meaning that the coefficient $\hat{\varphi}(x)$ is
a rational function).

Integrals such as (1) often appear as parametric representation of FIs. For
example, the Baikov representation of a FI with $\ell$ loops and $E$ external
legs in $d$ dimensions, the multivalued function $u(x)$ contains a single
factor $u(x)=\mathcal{B}(x)^{\gamma}$, where $\mathcal{B}$ is the Baikov
polynomial [25], and the exponent $\gamma=(d-\ell-E-1)/2$ . Hence in the
following we will refer to the integrals (1) as generalized Feynman Integrals
(GFI).

Linear equivalence relation between FIs:

$\displaystyle\mathcal{I}=\int_{\mathcal{C}}u\>\varphi\equiv\int_{\mathcal{C}}u\>(\varphi+\nabla_{\omega}\xi)\
,$ (2)

where we introduced the covariant derivative:
$\nabla_{\omega}:=\mathrm{d}+\omega\wedge$ and the $1$-form

$\displaystyle\omega:=\mathrm{d}\log{u}\ ,$ (3)

which will be very useful in the following. The equivalence relation (2)
follows from the Stokes theorem:
$0=\int_{\partial\mathcal{C}}u\>\xi=\int_{\mathcal{C}}u\>\nabla_{\omega}\xi\
,$ where $\nabla_{\omega}\xi:=\mathrm{d}\xi+\omega\wedge\xi$ is the covariant
derivative.

Fixing the contour of integration $\mathcal{C}$ allows us to interpret
relation (2) as an equivalence of integrands. Namely, we collect $n$-forms
$\varphi$ into equivalence classes
$\langle\varphi|:\varphi\sim\varphi+\nabla_{\omega}\xi$ generated by adding
covariant derivatives of $(n-1)$-forms. Their totality forms the twisted
cohomology group:

$\displaystyle\langle\varphi|\in\mathbb{H}^{n}_{\omega}:=\bigl{\\{}\text{$n$-forms
$\varphi$}\>|\>\nabla_{\omega}\varphi=0\bigr{\\}}\Big{/}\bigl{\\{}\nabla_{\omega}\xi\bigr{\\}}\
,$ (4)

which can be thought of as the space of linearly independent FIs (of a given
topology).

Analogously we can introduce the dual integrals $\mathcal{I}^{\vee}$, whose
definition mimics (1) up to $u\mapsto u^{-1}$ and
$\nabla_{\omega}\mapsto\nabla_{-\omega}$ . Elements of the dual twisted
cohomology group will be denoted by kets $|\psi\rangle$.

###  2.1 Counting the number of Master integrals

The framework of twisted cohomology unites several seemingly independent
methods for computation of the number of MIs $r$:

  1. 1.

Number of unreduced integrals produced by the Laporta algorithm [26].

  2. 2.

Number of critical points, i.e. solutions of $\mathrm{d}\log u(x)=0$ [27, 28].

  3. 3.

Number of independent integration contours $\mathcal{C}_{\lambda}$ [29, 30].

  4. 4.

Number of independent $n$-forms, i.e.
$\mathop{\mathrm{dim}}\bigl{(}\mathbb{H}^{n}_{\pm\omega}\bigr{)}$ [4, 8].

  5. 5.

Holonomic rank of GKZ system (volumes of certain polytopes) [1, 31].

###  2.2 Scalar product between Feynman integrals

The twisted cohomology theory allows us to view the set of FIs (of a given
topology) as a finite dimensional vector space. A set of MIs $\langle
e_{\lambda}|$ for $\lambda\in\\{1,\ldots,r\\}$ then forms a basis in that
space.

The dual FIs really form a dual vector space to FIs due to the existence of a
scalar product:

$\displaystyle\langle\varphi|\psi\rangle=\frac{1}{(2\pi\mathrm{i})^{n}}\int\iota(\varphi)\wedge\psi\
,$ (5)

called the intersection number. This scalar product allows to directly
decompose a given integral $\mathcal{I}$ in a basis of MIs
$\mathcal{J}_{\lambda}:=\int_{\mathcal{C}}u\>e_{\lambda}$, namely
$\mathcal{I}=\sum_{\lambda=1}^{r}c_{\lambda}\,\mathcal{J}_{\lambda}$ . Linear
algebra leads us to the master decomposition formula [4, 8]:

$\displaystyle\langle\varphi|=\sum_{\lambda=1}^{r}c_{\lambda}\,\langle
e_{\lambda}|\ ,\quad
c_{\lambda}=\sum_{\mu=1}^{r}\langle\varphi|h_{\mu}\rangle\bigl{(}C^{-1}\bigr{)}_{\mu\lambda}\
,$ (6) $\displaystyle C_{\lambda\mu}:=\langle e_{\lambda}|h_{\mu}\rangle\ ,$
(7)

for any choice of the dual basis $|h_{\mu}\rangle$. Therefore the intersection
numbers (5) completely determine the decomposition coefficients. Let’s see now
how they can be computed.

###  2.3 Univariate intersection numbers

In the $n=1$ case, intersection numbers (5) turn into a sum of residues [10,
11, 22]:

$\displaystyle\langle\varphi|\psi\rangle\equiv\frac{1}{2\pi\mathrm{i}}\int_{X}\left(\varphi-\sum_{p\in\mathbb{P}_{\omega}}\nabla_{\omega}\bigl{(}\theta_{p}(x,\bar{x})f_{p}\bigr{)}\right)\wedge\psi=\sum_{p\in\mathbb{P}_{\omega}}\mathop{\mathrm{Res}}_{x=p}\left[f_{p}\,\psi\right]\
,$ (8)

where

  * •

Integration goes over $X=\mathrm{\mathbb{C}P}^{1}$.

  * •

$\mathbb{P}_{\omega}:=\bigl{\\{}p\>\big{|}\>\text{poles of
$\omega$}\bigr{\\}}$, including the $\infty$ .

  * •

Terms with Heaviside $\theta$-functions regulate the integral with the help of
a local potential $f_{p}$, which satisfies
$\nabla_{\omega}f_{p}\equiv(\mathrm{d}+\omega\wedge)f_{p}=\varphi$ around the
pole $p$ . This differential equation can be solved via an Ansatz:
$f_{p}=f_{p,\,\mathrm{min}}(x-p)^{\mathrm{min}}+f_{p,\,\mathrm{min}+1}(x-p)^{\mathrm{min}+1}+\ldots+f_{p,\,\mathrm{max}}(x-p)^{\mathrm{max}}$
.

###  2.4 Multivariate intersection numbers

One strategy for dealing with the intersection numbers of multivariate FIs is
to apply the univariate procedure recursively one variable at a time [12, 18,
4, 7, 8]. Consider a 2 variable problem: given two 2-forms
$\varphi(x_{1},x_{2})$ and $\psi(x_{1},x_{2})$ we would like to compute
$\langle\varphi|\psi\rangle$ by first integrating out $x_{1}$ and then $x_{2}$
. To do that we pick a basis $\langle e_{\lambda}|$ and its dual
$|h_{\mu}\rangle$ for the internal $x_{1}$-integration and project $\varphi$,
$\psi$ onto them (omitting the summation signs):

$\displaystyle\langle\varphi|$ $\displaystyle=\langle
e_{\lambda}|\wedge\langle\varphi_{\lambda}|\ ,\quad$
$\displaystyle\langle\varphi_{\lambda}|=\langle\varphi|h_{\mu}\rangle\bigl{(}C^{-1}\bigr{)}_{\mu\lambda}\
,$ (9) $\displaystyle|\psi\rangle$
$\displaystyle=|h_{\mu}\rangle\wedge|\psi_{\mu}\rangle\ ,\quad$
$\displaystyle|\psi_{\mu}\rangle=\bigl{(}C^{-1}\bigr{)}_{\mu\lambda}\langle
e_{\lambda}|\psi\rangle\ ,$ (10)

The internal $x_{1}$-integration can be seen as the insertion of the identity
operator
$\mathbb{I}_{\mathrm{c}}=|h_{\mu}\rangle\bigl{(}C^{-1}\bigr{)}_{\mu\lambda}\langle
e_{\lambda}|\ ,$ which consequently allows us to write the remaining integral
in $x_{2}$ as a sum over residues:

$\displaystyle\langle\varphi|\psi\rangle=\langle\varphi\underbrace{|h_{\mu}\rangle\bigl{(}C^{-1}\bigr{)}_{\mu\lambda}\langle
e_{\lambda}|}_{\mathbb{I}_{\mathrm{c}}}\psi\rangle=\sum_{p\in\mathbb{P}_{P}}{\rm
Res}_{x_{2}=p}\left[f_{p,\lambda}\,C_{\lambda\mu}\,\psi_{\mu}\right]\ .$ (11)

Similar to Section 2.3, this formula requires the knowledge of a local vector
potential $f_{p,\lambda}$ near each pole $x_{2}=p$. The potential is fixed by
the following system of first order differential equations (omitting the $p$
subscript):

$\displaystyle\partial_{x_{2}}f_{\lambda}+f_{\mu}\>P_{\mu\lambda}=\varphi_{\lambda}\quad\text{near
$x_{2}=p$}\ .$ (12)

The differential equation matrix $P$ and it’s dual version $P^{\vee}$ are made
out of $x_{1}$-intersection numbers:

$\displaystyle
P_{\lambda\nu}:=\langle(\partial_{x_{2}}+\omega_{2})e_{\lambda}|h_{\mu}\rangle\bigl{(}C^{-1}\bigr{)}_{\mu\nu}\
,\quad P^{\vee}_{\mu\xi}:=\bigl{(}C^{-1}\bigr{)}_{\mu\lambda}\langle
e_{\lambda}|(\partial_{x_{2}}-\omega_{2})h_{\xi}\rangle\ ,$ (13)

so they can be computed using the univariate method of Section 2.3. The set
$\mathbb{P}_{P}$ in eq. (11) is defined as
$\mathbb{P}_{P}:=\bigl{\\{}p\>\big{|}\>\text{poles of $P$}\bigr{\\}}$ .

In practice, to compute the residue at, say, $x_{2}=0$ we solve for $\rho$ the
following system:

$\displaystyle\begin{cases}\left[x_{2}\,\partial_{x_{2}}+P(x_{2})\right]\vv{f}=\vv{\varphi}\\\\[5.0pt]
\rho=\mathop{\mathrm{Res}}_{x_{2}=0}\left[\vv{f}\,\cdot\,\vv{\psi}\right]\end{cases}\
,$ (14)

where we rescaled $\vv{\varphi}(x_{2})\mapsto 1/x_{2}\>\vv{\varphi}(x_{2})$
and $P(x_{2})\mapsto 1/x_{2}\>P(x_{2})$, and canceled the $C$ matrix in the
residue (11) against the $C^{-1}$ coming from eq. (10). The series expansion
of the system (14) is build from:

$\displaystyle P(x_{2})=\sum_{i\geq 0}x_{2}^{i}\>P_{i}\
,\quad\vv{\varphi}=\sum_{i\geq k}x_{2}^{i}\>\vv{\varphi}_{i}\
,\quad\vv{\psi}=\sum_{i\geq m}x_{2}^{i}\>\vv{\psi}_{i}\ ,$ (15)

for integer $k,m\in\mathbb{Z}$. Inserting an Ansatz
$\vv{f}=\sum_{i}x_{2}^{i}\>\vv{f}_{i}$ , and matching the powers of $x_{2}$
order by order, we obtain the linear system (here dots
${\color[rgb]{0.5,0.5,0.5}\cdot}$ denote zeros):

$\displaystyle\left[\begin{array}[]{c|ccccc|c}-1&\vv{\psi}_{1}&\vv{\psi}_{0}&\vv{\psi}_{-1}&\vv{\psi}_{-2}&\vv{\psi}_{-3}&\pagecolor{gr1}{\color[rgb]{0.5,0.5,0.5}\cdot}\\\
{\color[rgb]{0.5,0.5,0.5}\cdot}&P_{0}-2&{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}\cdot}&\vv{\varphi}_{-2}\\\
{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}P_{1}}&P_{0}-1&{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}\cdot}&\vv{\varphi}_{-1}\\\
{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}P_{2}}&{\color[rgb]{0.5,0.5,0.5}P_{1}}&P_{0}&{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}\cdot}&\vv{\varphi}_{0}\\\
{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}P_{3}}&{\color[rgb]{0.5,0.5,0.5}P_{2}}&{\color[rgb]{0.5,0.5,0.5}P_{1}}&P_{0}+1&{\color[rgb]{0.5,0.5,0.5}\cdot}&\vv{\varphi}_{1}\\\
{\color[rgb]{0.5,0.5,0.5}\cdot}&{\color[rgb]{0.5,0.5,0.5}P_{4}}&{\color[rgb]{0.5,0.5,0.5}P_{3}}&{\color[rgb]{0.5,0.5,0.5}P_{2}}&{\color[rgb]{0.5,0.5,0.5}P_{1}}&P_{0}+2&\vv{\varphi}_{2}\end{array}\right]\cdot\left[\begin{array}[]{c}\rho\\\
\vv{f}_{-2}\\\ \vdots\\\ \vv{f}_{2}\\\ -1\end{array}\right]=0\ .$ (27)

This equation has to be solved only for $\rho$ . Row reduction of this matrix
can be carried out only until the first row is filled with zeros except for
the element in the last column (highlighted with grey), which will contain the
needed residue. Other poles of eq. (11) are treated in the same manner and the
sum of their residues produces the intersection number
$\langle\varphi|\psi\rangle$ .

##  3 Decomposition via the secondary equation

As was observed in [1] (see also [31]), the twisted cohomology framework
provides another method for computation of the decomposition coefficients (7).
The first key idea is the so-called secondary equation [19, 8, 32], which is a
matrix differential equation satisfied by the intersection matrix $C$:

$\displaystyle\begin{cases}\partial_{z_{i}}\,\langle
e_{\lambda}|=\bigl{(}P_{i}\bigr{)}_{\lambda\nu}\,\langle e_{\nu}|\\\
\partial_{z_{i}}\,|h_{\mu}\rangle=|h_{\xi}\rangle\,\bigl{(}P^{\vee}_{i}\bigr{)}_{\xi\mu}\end{cases}\Longrightarrow\partial_{z_{i}}\,C=P_{i}\cdot
C+C\cdot\left(P_{i}^{\vee}\right)^{\mathrm{T}}\ ,$ (28)

where $z_{i}$ are some external kinematical variables. The other key step is
computation of the differential equation matrices $P$ and $P^{\mathrm{aux}}$
made available thanks to the connection of the twisted cohomology theory, the
GKZ formalism, and $D$-module theory. We assume that this step is completed
and refer the interested reader to [1, 31] for the full story. Once the
secondary equation (28) is written down, we employ the known algorithms for
finding rational solutions of such systems, e.g. the Maple package
IntegrableConnections [33].

Finally, to determine the decomposition coefficients (7) we repeat the above
procedure for an auxiliary basis
$e^{\mathrm{aux}}:=\\{e_{1},\ldots,e_{r-1},\varphi\\}$, i.e. we compute an
auxiliary $P^{\mathrm{aux}}$ and then $C^{\mathrm{aux}}$ . The FI
decomposition is then encoded in the following matrix product:

$\displaystyle\left[\begin{array}[]{c}e_{1}\\\ \vdots\\\ e_{r-1}\\\
\varphi\end{array}\right]=C^{\mathrm{aux}}\cdot
C^{-1}\left[\begin{array}[]{c}e_{1}\\\ \vdots\\\ e_{r-1}\\\
e_{r}\end{array}\right]\quad\Longrightarrow\quad C^{\mathrm{aux}}\cdot
C^{-1}=\left[\begin{array}[]{ccc|c}&&&0\\\ &{\mathrm{id}_{r-1}}&&\vdots\\\
&&&0\\\ \hline\cr c_{1}&\cdots&c_{r-1}&c_{r}\end{array}\right]\ ,$ (41)

where $\mathrm{id}_{r-1}$ denotes an identity matrix of size $(r-1)$, and the
decomposition coefficients $c_{\lambda}$ are collected in the last row
highlighted with grey.

###  3.1 A simple example

Let us briefly showcase how the secondary equation can produce the reduction
coefficients of a box diagram with a single dot $\varphi=\raisebox{-8.0pt}{
\axopicture{(20,20)(0,0)} \SetWidth{0.8} \Line(5,5)(15,5) \Line(15,5)(15,15)
\Line(15,15)(5,15) \Line(5,15)(5,5) \Line(0,0)(5,5) \Line(20,0)(15,5)
\Line(20,20)(15,15) \Line(0,20)(5,15) \Vertex(10,5){2} }$ in terms of the
basis $(e_{1},e_{2},e_{3})=\left(\raisebox{-8.0pt}{ \axopicture{(20,20)(0,0)}
\SetWidth{0.8} \ECirc(10,10){6} \Vertex(16,10){2} \Line(5,20)(10,16)
\Line(15,20)(10,16) \Line(5,0)(10,4) \Line(15,0)(10,4) },\raisebox{-8.0pt}{
\axopicture{(20,20)(0,0)} \SetWidth{0.8} \ECirc(10,10){6} \Vertex(10,16){2}
\Line(0,5)(4,10) \Line(0,15)(4,10) \Line(20,5)(16,10) \Line(20,15)(16,10)
\hskip 1.0pt },\raisebox{-8.0pt}{ \axopicture{(20,20)(0,0)} \SetWidth{0.4}
\Line(5,5)(15,5) \Line(15,5)(15,15) \Line(15,15)(5,15) \Line(5,15)(5,5)
\Line(0,0)(5,5) \Line(20,0)(15,5) \Line(20,20)(15,15) \Line(0,20)(5,15) }\
\right)$ . This topology has
$u=(x_{1}+x_{2}+x_{3}+x_{4}+x_{1}x_{3}+t/s\>x_{2}x_{4})^{\gamma}$ . Using the
algorithm of [1] and the Asir computer algebra system [34] we generate the
differential equation matrices:

$\displaystyle
P=\left[\begin{array}[]{ccc}-\frac{\epsilon\left(\delta^{2}(12z+11)+7\delta(z+1)+z+1\right)}{(3\delta+1)z(z+1)}&-\frac{\delta^{2}\epsilon}{(3\delta+1)(z+1)}&\frac{\delta^{2}\epsilon(\delta(z+2)+1)}{2(3\delta+1)z(z+1)(\delta\epsilon+1)}\\\
\frac{\delta^{2}\epsilon}{(3\delta+1)z\left(z+1\right)}&-\frac{\delta^{2}\epsilon}{(3\delta+1)(z+1)}&-\frac{\delta^{2}\epsilon(\delta+2\delta
z+z)}{2(3\delta+1)z(z+1)(\delta\epsilon+1)}\\\
-\frac{2(2\delta+1)\epsilon(\delta\epsilon+1)}{(3\delta+1)z(z+1)}&\frac{2(2\delta+1)\epsilon(\delta\epsilon+1)}{(3\delta+1)(z+1)}&-\frac{\epsilon\left(\delta^{2}(5z+7)+\delta(2z+5)+1\right)}{(3\delta+1)z(z+1)}\end{array}\right]\
,$ (45)

where $z=t/s$ is the ratio of the Mandelstam invariants, $\delta$ is an
additional regularization parameter which should be set $\delta\to 0$ at the
end of the computation, and $P^{\vee}=P\big{|}_{\epsilon\to-\epsilon}$ (see
[1, 31] for further details). The rational solution to the secondary equation
(28) looks like this:

$\displaystyle
C=\left[\begin{array}[]{ccc}-\frac{(2\delta+1)(4\delta+1)}{\delta}&\delta&-2(\delta\epsilon-1)\\\
\delta&-\frac{(2\delta+1)(4\delta+1)}{\delta}&-2(\delta\epsilon-1)\\\
2(\delta\epsilon+1)&2(\delta\epsilon+1)&-\frac{4\left(10\delta^{2}+6\delta+1\right)(\delta\epsilon-1)(\delta\epsilon+1)}{\delta^{3}}\end{array}\right]\
.$ (49)

We repeat the same procedure for the auxiliary $C^{\mathrm{aux}}$ and mount
everything into eq. (41) to produce:

$\displaystyle\raisebox{-8.0pt}{ \axopicture{(20,20)(0,0)} \SetWidth{0.8}
\Line(5,5)(15,5) \Line(15,5)(15,15) \Line(15,15)(5,15) \Line(5,15)(5,5)
\Line(0,0)(5,5) \Line(20,0)(15,5) \Line(20,20)(15,15) \Line(0,20)(5,15)
\Vertex(10,5){2}
}=-\frac{2\varepsilon(2\varepsilon+1)}{z(\varepsilon+1)}\cdot\raisebox{-8.0pt}{
\axopicture{(20,20)(0,0)} \SetWidth{0.8} \ECirc(10,10){6} \Vertex(16,10){2}
\Line(5,20)(10,16) \Line(15,20)(10,16) \Line(5,0)(10,4) \Line(15,0)(10,4)
}+0\cdot\raisebox{-8.0pt}{ \axopicture{(20,20)(0,0)} \SetWidth{0.8}
\ECirc(10,10){6} \Vertex(10,16){2} \Line(0,5)(4,10) \Line(0,15)(4,10)
\Line(20,5)(16,10) \Line(20,15)(16,10) \hskip 1.0pt
}+(2\varepsilon+1)\cdot\raisebox{-8.0pt}{ \axopicture{(20,20)(0,0)}
\SetWidth{0.4} \Line(5,5)(15,5) \Line(15,5)(15,15) \Line(15,15)(5,15)
\Line(5,15)(5,5) \Line(0,0)(5,5) \Line(20,0)(15,5) \Line(20,20)(15,15)
\Line(0,20)(5,15) }\ $

Therefore the secondary equation method allows us to decompose FIs in terms of
MIs via solving a first order matrix differential equation [1]!

##  4 Conclusion

We reviewed the connection between FIs and the twisted cohomology theory,
focusing on the algorithms for computation of the uni- and multivariate
intersection numbers, that is the scalar products between FIs.

Furthermore, following [1], we showed how the twisted cohomology together with
the theory of GKZ hypergeometric system provide a way to compute the IBP
reduction coefficients via essentially finding rational solutions to a system
of PDEs (28) called the secondary equation. In the future it would be
interesting to further develop this connection and apply it to other problems
and processes within pQFT.

Figures were made with AxoDraw2 [35].

## References

  * [1] V. Chestnov, F. Gasparotto, M.K. Mandal, P. Mastrolia, S.J. Matsubara-Heo, H.J. Munch et al., _Macaulay Matrix for Feynman Integrals: Linear Relations and Intersection Numbers_ , 2204.12983.
  * [2] F.V. Tkachov, _A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions_ , _Phys. Lett. B_ 100 (1981) 65.
  * [3] K.G. Chetyrkin and F.V. Tkachov, _Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops_ , _Nucl. Phys._ B192 (1981) 159.
  * [4] P. Mastrolia and S. Mizera, _Feynman Integrals and Intersection Theory_ , _JHEP_ 02 (2019) 139 [1810.03818].
  * [5] H. Frellesvig, F. Gasparotto, S. Laporta, M.K. Mandal, P. Mastrolia, L. Mattiazzi et al., _Decomposition of Feynman Integrals on the Maximal Cut by Intersection Numbers_ , _JHEP_ 05 (2019) 153 [1901.11510].
  * [6] S. Mizera, _Aspects of Scattering Amplitudes and Moduli Space Localization_ , 1906.02099.
  * [7] H. Frellesvig, F. Gasparotto, M.K. Mandal, P. Mastrolia, L. Mattiazzi and S. Mizera, _Vector Space of Feynman Integrals and Multivariate Intersection Numbers_ , _Phys. Rev. Lett._ 123 (2019) 201602 [1907.02000].
  * [8] H. Frellesvig, F. Gasparotto, S. Laporta, M.K. Mandal, P. Mastrolia, L. Mattiazzi et al., _Decomposition of Feynman Integrals by Multivariate Intersection Numbers_ , _JHEP_ 03 (2021) 027 [2008.04823].
  * [9] K. Matsumoto, _Quadratic Identities for Hypergeometric Series of Type $(k,l)$_, _Kyushu Journal of Mathematics_ 48 (1994) 335.
  * [10] K. Cho and K. Matsumoto, _Intersection theory for twisted cohomologies and twisted Riemann’s period relations I_ , _Nagoya Math. J._ 139 (1995) 67.
  * [11] K. Matsumoto, _Intersection numbers for logarithmic $k$-forms_, _Osaka J. Math._ 35 (1998) 873.
  * [12] K. Ohara, _Intersection numbers of twisted cohomology groups associated with Selberg-type integrals_ , 1998.
  * [13] K. Ohara, Y. Sugiki and N. Takayama, _Quadratic Relations for Generalized Hypergeometric Functions ${}_{p}F_{p-1}$_, _Funkcialaj Ekvacioj_ 46 (2003) 213.
  * [14] Y. Goto, _Twisted Cycles and Twisted Period Relations for Lauricella’s Hypergeometric Function $F_{C}$_, _International Journal of Mathematics_ 24 (2013) 1350094 [1308.5535].
  * [15] Y. Goto and K. Matsumoto, _The monodromy representation and twisted period relations for Appell’s hypergeometric function $F_{4}$_, _Nagoya Math. J._ 217 (2015) 61.
  * [16] Y. Goto, _Twisted period relations for Lauricella’s hypergeometric functions $F_{A}$_, _Osaka J. Math._ 52 (2015) 861.
  * [17] Y. Goto, _Intersection Numbers and Twisted Period Relations for the Generalized Hypergeometric Function ${}_{m+1}F_{m}$_, _Kyushu Journal of Mathematics_ 69 (2015) 203.
  * [18] S. Mizera, _Scattering Amplitudes from Intersection Theory_ , _Phys. Rev. Lett._ 120 (2018) 141602 [1711.00469].
  * [19] S.-J. Matsubara-Heo and N. Takayama, _An algorithm of computing cohomology intersection number of hypergeometric integrals_ , _Nagoya Mathematical Journal_ (2019) 1 [1904.01253].
  * [20] P. Mastrolia, _From Diagrammar to Diagrammalgebra_ , _PoS_ MA2019 (2022) 015.
  * [21] S. Mizera, _Status of Intersection Theory and Feynman Integrals_ , _PoS_ MA2019 (2019) 016 [2002.10476].
  * [22] H.A. Frellesvig and L. Mattiazzi, _On the Application of Intersection Theory to Feynman Integrals: the univariate case_ , _PoS_ MA2019 (2022) 017 [2102.01576].
  * [23] M.K. Mandal and F. Gasparotto, _On the Application of Intersection Theory to Feynman Integrals: the multivariate case_ , _PoS_ MA2019 (2022) 019.
  * [24] S. Weinzierl, _Feynman Integrals_ , 2201.03593.
  * [25] P.A. Baikov, _Explicit solutions of the multiloop integral recurrence relations and its application_ , _Nucl. Instrum. Meth. A_ 389 (1997) 347 [hep-ph/9611449].
  * [26] S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys. A_ 15 (2000) 5087 [hep-ph/0102033].
  * [27] P.A. Baikov, _A Practical criterion of irreducibility of multi-loop Feynman integrals_ , _Phys. Lett. B_ 634 (2006) 325 [hep-ph/0507053].
  * [28] R.N. Lee and A.A. Pomeransky, _Critical points and number of master integrals_ , _JHEP_ 11 (2013) 165 [1308.6676].
  * [29] J. Bosma, M. Sogaard and Y. Zhang, _Maximal Cuts in Arbitrary Dimension_ , _JHEP_ 08 (2017) 051 [1704.04255].
  * [30] A. Primo and L. Tancredi, _Maximal cuts and differential equations for Feynman integrals. An application to the three-loop massive banana graph_ , _Nucl. Phys. B_ 921 (2017) 316 [1704.05465].
  * [31] H.J. Munch, _Feynman Integral Relations from GKZ Hypergeometric Systems_ , _PoS_ LL2022 (2022) .
  * [32] S. Weinzierl, _On the computation of intersection numbers for twisted cocycles_ , _J. Math. Phys._ 62 (2021) 072301 [2002.01930].
  * [33] M.A. Barkatou, T. Cluzeau, C. El Bacha and J.-A. Weil, _Computing closed form solutions of integrable connections_ , in _Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’12, p. 43–50, 2012, DOI.
  * [34] “Risa/Asir, OpenXM project.” http://www.math.kobe-u.ac.jp/Asir/.
  * [35] J.C. Collins and J.A.M. Vermaseren, _Axodraw Version 2_ , 1606.01177.


# Generating functions in symplectic and contact geometry

A translated point of a contactomorphism $\phi$ on a contact manifold with
contact form $\alpha$ is a point $p$ where $\alpha$ is preserved under $\phi$
and whose image under $\phi$ lies in the same Reeb trajectory. They were
introduced as a contact analogon for fixed points of Hamiltonian
diffeomorphisms by Sheila Sandon in [San12] and can be understood as a special
case of leafwise fixed points. A contact version of the non-degenerate Arnol’d
conjecture on spheres was established in [San13] using a generating function
approach. It turns out that Sandon’s proof only works under the assumption
that there exists a generating function whose sublevel set at zero has
nontrivial homology. This thesis proves the result under this additional
assumption and fills gaps in other parts of Sandon’s argument.

figure . definition]Proposition definition]Corollary definition]Theorem
definition]Lemma definition]Def. & Prop definition]Conjecture 3.6cm3.7cm*
3.9cm2.7cm* dash subsection

Generating functions in
symplectic and contact geometry

Master’s thesis by
Aaron Gootjes-Dreesbach

submitted to |  Department of Mathematics
---|---
|  Universiteit Utrecht
supervised by |  Dr. Fabian Ziltener
second reader |  Prof. Dr. Marius Crainic

September 2020

Updated version from .

Abstract

Contents

*


# Generation of 10-m-lengthscale plasma columns by resonant and off-resonant
laser pulses

G. Demeter  J. T. Moody111Now at: Lawrence Livermore National Laboratory,
California, USA M. Á. Kedves  F. Batsch  M. Bergamaschi  V. Fedosseev  E.
Granados  P. Muggli  H. Panuganti  G. Zevi Della Porta  Wigner Research Centre
for Physics, Budapest, Hungary  Max Planck Institute for Physics, Munich,
Germany  CERN, Geneva, Switzerland

###### Abstract

Creating extended, highly homogeneous plasma columns like that required by
plasma wakefield accelerators can be a challenge. We study the propagation of
ultra-short, terawatt power ionizing laser pulses in a 10-meter-long rubidium
vapor and the plasma columns they create. We perform experiments and numerical
simulations for pulses with 780 nm central wavelength, which is resonant with
the D2 transition from the ground state of rubidium atoms, as well as for
pulses with 810 nm central wavelength, some distance from resonances. We
measure transmitted energy and transverse width of the pulse and use schlieren
imaging to probe the plasma column in the vapor close to the end of the vapor
source. We find, that resonant pulses are more confined in a transverse
direction by the interaction than off-resonant pulses are and that the plasma
columns they create are more sharply bounded. Off-resonant pulses leave a
wider layer of partially ionized atoms and thus lose more energy per unit
propagation distance. Using experimental data, we estimate the energy required
to generate a 20-meter-long plasma column and conclude that resonant pulses
are much more suitable for creating a long, homogeneous plasma.

###### keywords:

High-power pulse propagation , Resonant nonlinear interaction , Schlieren
imaging , Plasma wakefield acceleration

##  1 Introduction

The propagation of high-power, ionizing laser pulses in gases has been studied
extensively for decades, under a wide range of conditions [1, 2, 3]. Phenomena
such as self-focusing, filamentation or super-continuum generation were
investigated and numerous applications for these phenomena devised (remote
sensing, nonlinear spectroscopy, lightning protection, etc.). One particular
application is the creation of plasma columns to be used in wakefield particle
acceleration. Plasma wakefield accelerators are capable of accelerating
electrons (or positrons) in the strong electric fields sustained by plasma
waves, [4, 5, 6, 7, 8]. Since the accelerating gradients can be up to three
orders of magnitude larger than in conventional particle accelerators,
wakefield acceleration may be a replacement for established accelerator
technology in compact, cost-effective particle accelerators for science and
commercial applications. A multitude of scientific and technological
challenges need to be tackled (only one of which is the creation of the
necessary plasma column), but wakefield accelerators and their applications
are advanced constantly by numerous research groups worldwide [9, 10, 11, 12].
One notable example of such a wakefield accelerator project is the Advanced
Wakefield Experiment (AWAKE) at CERN (European Laboratory for Particle Physics
in Geneva, Switzerland), the first proton driven wakefield accelerator [13,
14].

At the heart of the AWAKE device is a 10-meter-long plasma column with finely
engineered plasma density, which is essential for accelerator operation. High
energy proton bunches from the Super Proton Synchrotron facility at CERN
interact with this plasma to create large amplitude wakefields, which in turn
can accelerate electron bunches well into the GeV domain. This plasma column
is created by starting from a vapor of rubidium atoms with precisely tuned
density distribution along the vapor source axis and propagating a terawatt
(TW) power laser pulse along the vapor to achieve exactly one-electron
ionization of atoms with a probability very close to one. Precisely tuned
vapor density thus yields a finely tuned plasma density. Achieving complete
one-electron ionization is facilitated by the fact that the 780 nm wavelength
laser pulse is resonant with the
$\mathrm{5s^{2}S_{1/2}}\rightarrow\mathrm{5p^{2}P_{3/2}}$ transition of the
valence electron from the atomic ground state (the D2 line) and then further
from $\mathrm{5p^{2}P_{3/2}}$ to $\mathrm{5d^{2}D_{5/2}}$,
$\mathrm{5d^{2}D_{3/2}}$ states. These single-photon resonances also have a
major impact on the propagation of the ionizing laser pulse – a question
studied only very recently in the context of ultra-short pulse propagation
[15]. It has been suggested [16], that it is predominantly the single-photon
resonances that give rise to a strong but saturable nonlinearity which can be
very advantageous for the propagation of the ionizing pulse. In particular,
theoretical indication was that due to resonant self-focusing of the pulse by
the vapor, pulse energy in the tail is channeled very effectively within the
plasma created by the leading edge of the pulse. The plasma column therefore
becomes longer and more sharply bounded when a 780 nm ionizing pulse is used.

It is important to note that the laser pulse in this scheme is not intense
enough to drive wakefields in the plasma (such as in laser wakefield
acceleration). Field strengths are orders of magnitude smaller and therefore
associated nonlinear propagation effects (such as relativistic self-focusing)
are absent. The scheme is related to plasma column formation for laser
wakefield acceleration with the help of pre-formed plasma waveguides [17, 18].
However, as plasma densities are orders of magnitude smaller, inhomogeneities
of plasma dispersion itself play little role during propagation.

Here we present an experimental investigation of plasma column generation by
resonant and off-resonant TW laser pulses. We compare the propagation
properties of pulses with 780 nm central wavelength (the rubidium D2 line) and
pulses with 810 nm central wavelength in the 10-meter-long vapor source at the
CERN AWAKE site. We also measure the properties of the created plasma column
close to the downstream end of the vapor source using schlieren imaging. Using
measurement data, an extension of the plasma column creation with the same
laser apparatus to 20 meter length is considered. We generalize the theory
derived for resonant pulse propagation [15] to treat both resonant and off-
resonant pulses with the same equations and we perform numerical simulations
to study the pulse propagation process. Comparing measurement results with
calculations we show that the predictions of our theory are qualitatively
correct for both resonant and off-resonant propagation. However, some
quantitative discrepancies between simulation and experiment remain in certain
respects. The primary application of our results is in wakefield accelerator
design, but they are also interesting for any application associated with the
propagation of high-power laser pulses such as the creation of long plasma
channels for lightning protection [19] or remote sensing applications where
resonances play an important role.

##  2 Experiment

###  2.1 Laser propagation experiment apparatus

Experiments were performed at the AWAKE site at CERN, with the rubidium vapor
source of the wakefield accelerator device [13, 20, 21], a 10 m long, 4 cm
diameter, temperature controlled steel tube that contains the Rb vapor. A
schematic drawing of the experimental setup can be seen on Fig. 1. A TW class
Ti:Sa laser system supplied $\sim$120 fs duration, $\sim$150 mJ energy pulses
for the ionization of the vapor. A mismatched, $\sim$40 m effective focal
length telescope was used to focus the pulses into the vapor source through a
10-mm-diameter aperture. The waist diameter was $w\approx 1.8$ mm full width
at half maximum (FWHM), the waist location was tuned near the center of the
vapor source. In half of the measurements, spectral shaping methods were
implemented to confine the spectrum of a wide bandwidth Ti:Sa oscillator to a
region around 780 nm wavelength, precisely the wavelength of the Rb D2
resonance line, the $\mathrm{5s^{2}S_{1/2}}\rightarrow\mathrm{5p^{2}P_{3/2}}$
transition from the atomic ground state (Fig. 2). The spectrum also had
significant intensity at the 795 nm $D_{1}$ resonance line
($\mathrm{5s^{2}S}_{1/2}\rightarrow\mathrm{5p^{2}P}_{1/2}$ transition) and the
776 nm transitions to higher lying excited states
($\mathrm{5p^{2}P}_{3/2}\rightarrow\mathrm{5d^{2}D}_{3/2}$ and
$\mathrm{5d^{2}D}_{5/2}$), similarly to our previous experiments [16]. In the
other half of the measurements, spectral shaping was used to obtain laser
pulses with the central wavelength shifted away from these resonances. The
resulting amplified pulse spectrum had a central wavelength of 810 nm and
virtually no power at the Rb resonance wavelengths (Fig. 2). These off-
resonant pulses had about 20% less maximum input energy than the 780 nm
central wavelength resonant pulses.

Precise temperature control of the vapor source reservoir and walls made it
possible to create a regulated, homogeneous, constant density Rb vapor
($\delta\rho/\rho<0.5\%$) in a $\mathcal{N}=10^{14}-10^{15}$ cm-3 range. We
measured rubidium density at the upstream end of the vapor source using white-
light interferometry [22, 23, 24]. The laser pulse energy was regulated by a
half-waveplate and two thin-film polarizers at Brewster angle between the last
amplifier and the compressor. A “virtual” laser line was set up using the
transmission from one of the transport mirrors in the laser line upstream of
the vapor source. This line had cameras that recorded the laser pulse
transverse energy distribution at propagation distances corresponding to the
entrance (VLC1), center (VLC2) and exit (VLC3) of the vapor source (Fig. 1).
These images represent propagation of the laser pulse in vacuum, for
comparison with its propagation through the vapor source and for monitoring
the focusing. An energy meter was also placed in the virtual line to measure
pulse input energy $E_{in}$, calibrated using a direct energy meter when the
vapor source vacuum system was open.

Figure 1: Sketch of the experimental setup. Pulses from a TW laser system
(Ti:Sa) propagate along the rubidium vapor source (Rb). About 1% of the pulse
energy is deflected to the virtual laser line to be monitored by an energy
meter ($E_{in}$) and cameras (VLC1-3). About 0.5% of the transmitted pulse is
reflected off a wedge to an energy meter ($E_{out}$) and through an imaging
telescope (T) to a camera (Pickoff). A transverse probe beam close to the
downstream end of the vapor source is used for schlieren imaging on a gated
camera (GC).

Downstream of the vapor source, a wedge before the beam dump diverted $\sim
0.5$% of the transmitted laser pulse to the output energy meter $E_{out}$ and
an imaging system that created an image of the vapor source output aperture on
the pickoff camera. This was used to record the ionizing pulse transverse
energy profile after propagating through the vapor. We calibrated $E_{out}$
readings to $E_{in}$ values by a series of measurements performed with only
residual rubidium vapor in the chamber ($\mathcal{N}\ll 0.5\cdot 10^{14}$
cm-3), which is estimated to absorb energy from the pulse well below the
$\sim$ mJ noise floor of the output energy measurement and affect negligibly
the pulse energy distribution. We used the same measurements to scale the size
of the pickoff images to the known size of the virtual exit camera image.
Various filters were used on each of the virtual laser line cameras and the
pickoff camera to prevent image saturation. The supplemental material of
reference [16] contains a detailed drawing of the experimental setup for the
energy and transmitted pulse measurement, the laser virtual line and a
description for the calibration procedure for the output energy and
transmitted pulse transverse profile.

Figure 2: Measured spectrum of the resonant (red) and off-resonant (blue)
laser pulses. Single-photon resonances to the $\mathrm{5p^{2}P_{1/2}}$ (795
nm) and $\mathrm{5p^{2}P_{3/2}}$ (780 nm) first excited states are marked by
black lines, as well as the resonance from $\mathrm{5p^{2}P_{3/2}}$ to higher
lying excited states $\mathrm{5d^{2}D_{3/2}}$ and $\mathrm{5d^{2}D_{5/2}}$
(776 nm - see also Fig. 11).

In order to examine the difference between the interaction of resonant and
off-resonant laser pulses with the Rb vapor, we performed propagation
measurements for several vapor densities, switching the laser spectrum each
time to take data with both resonant and off-resonant ionizing laser pulses.
Measurements were done for $\mathcal{N}\approx 2\cdot
10^{14}\mathrm{~{}cm^{-3}}$, $\mathcal{N}\approx 5\cdot
10^{14}\mathrm{~{}cm^{-3}}$, $\mathcal{N}\approx 7\cdot
10^{14}\mathrm{~{}cm^{-3}}$ and $\mathcal{N}\approx 1\cdot
10^{15}\mathrm{~{}cm^{-3}}$ vapor densities, the precise value differing by
less than $\pm 2\%$ for the corresponding resonant pulse / off-resonant pulse
measurements. For $\mathcal{N}\approx 5\cdot 10^{14}\mathrm{~{}cm^{-3}}$ and
$\mathcal{N}\approx 7\cdot 10^{14}\mathrm{~{}cm^{-3}}$ vapor densities, we
also measured the properties of the plasma column close to the downstream end
of the vapor source using schlieren imaging together with the transmitted
pulse measurements.

###  2.2 Schlieren imaging of the plasma column

Schlieren imaging is a very sensitive method to measure refraction index
changes in transparent media, used predominantly in aeronautics and fluid
mechanics [25], but also employed regularly to investigate laser induced
plasma [26, 27, 28, 29, 30]. Recently it was tested to probe atomic excitation
in rubidium vapor [31] and to measure plasma column properties in ionized
rubidium [32]. In our setup, the $\lambda_{p}=780.311$ nm (in vacuum)
continuous wave diode laser probe, tuned close to the precise value of the
$\lambda_{D2}=780.241$ nm (in vacuum) D2 resonance line, crossed the vapor
source in a transverse direction through a pair of sapphire view ports, as
seen on Fig. 3. The vapor refractive index contribution is thus $\delta
n=10^{-4}-10^{-3}$, due to anomalous dispersion by ground state atoms. This is
large enough for detection despite the low density of vapor compared to
standard atmospheric densities and it is 3-4 orders of magnitude larger than
the refractive index contribution due to plasma dispersion at the same
density. When the ionizing laser creates plasma, the population of the atomic
ground state is reduced significantly, so the refractive index changes
locally. A spatially dependent phase shift is imprinted upon the probe beam
during its transit.

Figure 3: a) Sketch of the schlieren imaging setup. A probe beam (PB)
traverses the vapor source (C) cross section through a pair of sapphire view
ports (W) along the $z$ direction, sampling the refractive index of the vapor
and the plasma column (P). With two lenses (L1, L2) in a $4f$ setup and a mask
(M) between them, a schlieren image is created on the gated camera (GC). b) A
narrow region of the gated camera image is extracted and averaged to obtain
the signal c), which is then frequency filtered before being evaluated d).

We placed two 75 cm focal length lenses in a $4f$ setup [33] after the vapor
cell and circular mask with a diameter of $D_{m}=1.5$ mm at the back focal
plane of the first lens. We used a gated, image-intensified camera (Andor
iStar DH334T-18F-73) to retrieve images of the probe light, triggered 100 ns
after the ionizing pulse and timed to collect light for a duration of 100 ns.
With this timing, atoms excited, but not ionized by the laser will return to
the ground state by spontaneous decay, while plasma recombination (1-10 $\mu$s
timescale at these densities) will not alter the plasma density considerably
yet. Thus probing the absence of ground state atoms yields information on the
plasma density.

A typical schlieren image captured by the measurement can be seen on Fig. 3
b). Given that the properties of the plasma column are constant on the mm
scale along the ionizing pulse propagation direction $x$, the $w=3.1$ mm FWHM
diameter probe beam samples a $y-z$ dependent vapor density, the $x$
dependence on the image is due only to the probe beam intensity variation.
Therefore it is convenient to take a rectangular region from the image, narrow
in the $x$ direction around the center of the probe beam and average along $x$
to improve the signal-to-noise ratio. The 1D $s(y)$ curve obtained is the
schlieren signal that we use to analyze the plasma column (shown in Fig. 3 c)
).

###  2.3 Inferring plasma column properties

The parameters of the schlieren probe beam were determined in a series of
measurements with the mask removed and only residual vapor density in the
vapor source. Given these parameters and the precise data of the anomalous
dispersion of the vapor [34], the schlieren signal $s(y)$ can also be
calculated theoretically using standard formulas of Fourier optics [33] for
any given plasma density distribution $N_{p}(y,z)$. (Note however, that at the
vapor densities considered, the homogeneous lineshape function in [34] must be
augmented by a collision broadening term [35], the magnitude of which contains
a constant known experimentally to much lower accuracy than the spontaneous
decay rate.) To obtain information on the extent of the plasma column from
$s(y)$, we start by assuming some sensible profile for the plasma density and
calculating the theoretical schlieren signal. According to theory [15, 16] and
experiment [36], when the ionizing pulse is powerful enough, the ionization
fraction at the center of the column is very close to one, i.e. plasma density
is saturated at the initial vapor density. Thus we assume an axisymmetric
plasma density of the form:

${N}_{p}(y,z)=\left\\{\begin{aligned}
&\mathcal{N}P_{max},&\mathrm{~{}if~{}}r\leq r_{0},\\\
&\mathcal{N}P_{max}\exp\left(-\frac{(r-r_{0})^{2}}{t_{0}^{2}}\right),&\mathrm{~{}if~{}}r>r_{0}.\end{aligned}\right.$
(1)

Here $\mathcal{N}$ is the rubidium vapor density, $r$ is the distance from the
column center, which is located at coordinates $(y_{0},z_{0}=0)$ in the $y-z$
plane. (A nonzero value of $z_{0}$ does not change the schlieren signal
because probe light phase modulation and absorption arise as dielectric
parameter integrals along $z$.) $r_{0}$ is the radius of the plasma column
core, $P_{max}$ is the maximum ionization fraction at the center and $t_{0}$
is the sheath layer width parameter, a value that characterizes the width of
the transition region between the core and the completely unionized vapor of
neutral atoms. When the vapor is not ionized completely at the center
($P_{max}<0$), we expect $r_{0}\approx 0$, and a clear maximum of ionization
fraction at the center. A substantial value of $r_{0}$ is only compatible with
$P_{max}\approx 1$, the saturation of ionization in the core. The reason is
that an extended region of constant ionization fraction can be realized either
if we realize a sizable region in the transverse plane where the time-
dependent laser pulse intensity is the same, or if the fluence is simply high
enough for the ionization fraction to saturate to 1. The former is very hard
to imagine in a system with a propagating pulse (self-focusing, diffraction).

The choice of the plasma profile Eq. 1 (in particular the Gaussian decay
function outside the core) is motivated by the fact that in the limit of less
intense fields (multiphoton ionization) and pulse profile undistorted by
nonlinear self-focusing (at the entrance of the vapor cell), a Gaussian beam
profile will give rise to a Gaussian plasma profile. For the general case,
(pulse profile already distorted by the nonlinear interaction during
propagation) the sheath layer width parameter $t_{0}$ can simply be regarded
as a parameter of a function fit. It will depend on the transverse fluence
distribution of the propagating laser pulse and ultimately it helps
characterize the strength of the self-focusing effect of the vapor. Note that
this width is not related to the plasma Debye length, as the ionization takes
place on the 100 fs timescale, much shorter than the timescale for any plasma
dynamical phenomena.

Calculating the schlieren signal $s(y)$ for plasma columns described by Eq. 1
with a range of sensible parameter values, one can verify that the plasma
gives rise to a double peaked structure in $s(y)$, similar to the experimental
signal (Fig. 3). The two peaks are due to the probe light phase modulation
varying in space most near the top and bottom edges of the plasma column. It
is convenient to frequency-filter $s(y)$ numerically with a low pass filter to
remove spatial frequencies $f_{y}\geq D_{mask}/(\lambda_{p}l)$, as the
interference of light passing above and below the mask edges distorts the
peaks somewhat ($l=75$ cm is the distance between the mask and the second
lens). We can compute the locations $y_{1},y_{2}$ of the two largest peaks of
the filtered schlieren signal $\mathcal{F}\left(s(y)\right)$, as well as the
peak widths $\sigma_{1},\sigma_{2}$ and peak heights $A_{1},A_{2}$ using (a
slightly tweaked version of) the find_peaks function of SciPy [37]. It is then
possible to verify that the core center location $y_{0}$ can be determined by
$y_{0}=(y_{1}+y_{2})/2$ with very good accuracy. Furthermore, defining the
peak distance $\Delta$ and normalized peak width $W$:

$\displaystyle\Delta$ $\displaystyle=$ $\displaystyle|y_{2}-y_{1}|$ (2)
$\displaystyle W$ $\displaystyle=$
$\displaystyle\frac{A_{1}\sigma_{1}+A_{2}\sigma_{2}}{A_{1}+A_{2}}$ (3)

we observe that these quantities are, to a very good approximation, linear
functions of $r_{0}$ for fixed $t_{0}$ and vice versa, while they do not
depend on $y_{0}$ at all. Therefore we write their functional dependence in
the following form:

$\displaystyle\Delta$ $\displaystyle=M_{12}r_{0}t_{0}+M_{1}r_{0}+M_{2}t_{0}+B$
$\displaystyle W$ $\displaystyle=Q_{12}r_{0}t_{0}+Q_{1}r_{0}+Q_{2}t_{0}+P$ (4)

and use a set of signals calculated with varying $r_{0}$, $t_{0}$ and $y_{0}$
to determine, using a fitting procedure, the set of constants
$\\{M_{1},M_{2},M_{12},B,Q_{1},Q_{2},Q_{12},P\\}$ from $\Delta(r_{0},t_{0})$
and $W(r_{0},t_{0})$. Once known, we can evaluate any schlieren signal we
obtain from the experiment by spatial frequency filtering to get rid of
interference (and high-frequency noise components) and using peak finding to
determine $y_{0}$, $\Delta$ and $W$. We can then associate with the signal the
plasma core radius $r_{0}$ and the sheath layer width $t_{0}$ obtained by
inverting the relations Eqs. 4. In a sense, we can regard $r_{0}$ and $t_{0}$
as the parameters of a function fit of the form Eq. 1 on the experimental
plasma distribution.

For some parameter combinations, the two-peaked structure may be absent, so we
cannot associate $r_{0},t_{0}$ values with $s(y)$ using the procedure — in
this case the image cannot be evaluated. This may happen if sufficiently large
values of $y_{0}$ and/or $r_{0}$ and/or $t_{0}$ combine such that one or both
of the plasma edges lie far to one side where the probe beam is already too
weak. This may also happen if $r_{0}$ and $t_{0}$ are both small and the two
peaks are not separated. With the present measurement, $r_{0}+t_{0}\gtrsim
0.3$ mm is required for a reliable separation of the peaks. While frequency
filtering the schlieren signal to mitigate interference effects increases the
accuracy of the evaluation, it carries a price. Very sharp plasma boundaries
($t_{0}\leq 0.1$ mm) give rise to narrow peaks in the schlieren image and are
distorted by the filtering we use. This effectively sets the lower limit on
the sheath width we can reliably evaluate. Any $t_{0}$ below this limit will
be measured as $t_{0}\approx 0.1$ mm.

Note that $\Delta$ and $W$ do not depend on the overall magnitude of the
signal. This is convenient because vapor absorption is not known to a high
accuracy due to collision broadening and also because probe laser power was
not monitored continuously. The evaluation process just described yields no
information on $P_{max}$ – but when there is a measurable plasma column core
($r_{0}>0$) we can safely assume $P_{max}=1$. We further note that the fit
coefficients (with given beam parameters) depend somewhat on the vapor
density, so for the evaluation of any measurement, the corresponding set of
theoretical samples must be computed and the fitting parameters determined.

Estimating the plasma parameters using machine learning methods has also been
tested previously [38]. Deep neural networks were trained using a large number
of calculated signals to recover the underlying parameters and they proved
more accurate than the fitting procedure in this paper for the calculated
signals. However, for the actual experimental data, their predictions
exhibited large fluctuations at times, most probably because real plasma
columns are not axially symmetric and may have slightly different sheath
thickness at the top and the bottom edge. The present evaluation method proves
more robust with respect to this circumstance.

##  3 Resonant vs. off-resonant pulse propagation

###  3.1 Experimental observations

As a representative example of the measurement results for the transmitted
pulse properties we obtained, Fig. 4 shows the width of the transmitted pulse
together with its pulse $E_{out}$ as a function of $E_{in}$, measured for
$\mathcal{N}=7\cdot 10^{14}\mathrm{~{}cm}^{-3}$ vapor density, resonant ($780$
nm) pulses. To characterize the beam width, we use the D4$\sigma$ width (i.e.
the second moment width) of the fluence profile $\mathcal{F}(x,y)$ defined as:

$D4\sigma=4\sqrt{\frac{\int\left(\mathcal{F}(x,y)(x-\bar{x})^{2}+\mathcal{F}(x,y)(y-\bar{y})^{2}\right)dxdy}{\int\mathcal{F}(x,y)dxdy}}$
(5)

where

$\bar{x}=\frac{\int
x\mathcal{F}(x,y)dxdy}{\int\mathcal{F}(x,y)dxdy}\mathrm{~{}and~{}}\bar{y}=\frac{\int
y\mathcal{F}(x,y)dxdy}{\int\mathcal{F}(x,y)dxdy}.$ (6)

$x$ and $y$ are the coordinates in the camera plane in these formulas,
$\bar{x},\bar{y}$ are coordinates of the geometric center and clearly
$\int\mathcal{F}(x,y)dxdy=E_{out}$. The figure depicts “raw” measurement data,
each marker corresponds to a single measurement. Insets show transmitted pulse
camera pictures for a few single representative measurements. Several regimes
are visible on the plots, as discussed in [16]. For small $E_{in}$ in the sub-
threshold domain (marked “ST” on Fig. 4), $E_{out}$ is below the noise-floor
of the measurement and the image on the pickoff camera is very broad (see
inset a) ). Then follows the sharp breakthrough transition (marked “B” on Fig.
4), where transmitted pulse width drops quickly as it develops a sharp,
narrow, high-fluence feature and $E_{out}$ starts to increase. Just above this
is the confined beam domain (“CB” on Fig. 4), where the pulse width slowly
approaches a minimum, while $E_{out}$ increases (insets b) and c) ). In the
final region named the asymptotic transparency domain (“AT” on Fig. 4), both
$E_{out}$ and pulse width increase steadily (inset d) ). This region is
associated with the saturation of the optical nonlinearity of the medium due
to complete one-electron ionization in the plasma column core ($P_{max}=1$).

Figure 4: Measurement results for transmitted pulse D4$\sigma$ width (right
$y$-axis, red crosses) and energy $E_{out}$ (left $y$-axis, blue plus symbols)
as a function of pulse input energy for resonant 780 nm pulses and
$\mathcal{N}=7\cdot 10^{14}\mathrm{~{}cm}^{-3}$ vapor density. Dashed vertical
lines mark approximate domain boundaries, insets show representative
transmitted pulse shapes: a) pulse in the sub-threshold (ST) domain, b)
narrow-width pulses in the breakthrough (B) and c) confined beam (CB) domains
and finally d) widening pulse of the asymptotic transparency (AT) domain.

The main results of the propagation experiments can be seen on Figs. 5 where
$E_{out}$ and transmitted pulse D4$\sigma$ width can be seen for all four
vapor densities studied and both ionizing pulse central wavelengths. Binned
data averages are plotted with error bars marking the standard error of the
mean. The figures show that off-resonant ionizing pulses (blue lines) behave
similarly to resonant pulses (red lines) in general. However, $E_{out}$ is
smaller for any given $E_{in}$. The breakthrough transition requires higher
$E_{in}$ and the transmitted pulse is always wider for off-resonant pulses. As
the vapor density increases, the breakthrough transition shifts to higher
$E_{in}$ for both wavelengths. Apart from the lowest density measurements, the
off-resonant pulse also acquires a minimum width after the breakthrough. The
larger $E_{out}$ and smaller D4$\sigma$ width combine to give rise to a
substantially larger peak fluence in the resonant case than in the off-
resonant case. Overall, Figs. 5 b), d), f) and h) show that the nonlinear
self-focusing effect of the vapor is stronger on the resonant pulses than on
the off-resonant pulses, in accordance with theoretical predictions.

Figure 5: Transmitted pulse energy $E_{out}$ (left panels) and D4$\sigma$
width (right panels) as a function of pulse energy for all four vapor
densities and both ionizing pulse wavelengths. Markers show binned data
averages, error bars depict the standard error of the mean. Red symbols/lines
mark 780 nm resonant, blue lines/symbols mark 810 nm off-resonant
measurements. a), b) $\mathcal{N}=2\cdot 10^{14}\mathrm{~{}cm}^{-3}$; c), d)
$\mathcal{N}=5\cdot 10^{14}\mathrm{~{}cm}^{-3}$; e), f) $\mathcal{N}=7\cdot
10^{14}\mathrm{~{}cm}^{-3}$; g), h) $\mathcal{N}=1\cdot
10^{15}\mathrm{~{}cm}^{-3}$.

We used schlieren imaging to measure plasma column dimensions for the two
intermediate vapor densities together with the propagation measurement. Figure
6 shows the plasma core radius $r_{0}$ we obtained, as a function of $E_{in}$.
Transmitted pulse D4$\sigma$ width is also shown on the figure for reference.
It is visible that a measureable $r_{0}$ appears at the breakthrough
transition where the transmitted pulse width drops drastically. For resonant
pulses, following an initial sharp increase of $r_{0}$ there is a visible
“shoulder” of near constant plasma radius, approximately corresponding to the
confined beam domain ($E_{in}=32-50\mathrm{~{}mJ}$ and
$E_{in}=45-72\mathrm{~{}mJ}$ for the two vapor densities depicted). There is
no such feature visible for off-resonant pulses. The plasma column reaches the
downstream end of the 10 meter vapor source at higher $E_{in}$ for off-
resonant pulses (as the breakthrough requires higher energies), there is a
roughly $r_{0}=0.5-0.6$ mm plasma radius with resonant pulses already when the
first detectable plasma appears for the off-resonance case. However, by the
end of the input energy range both resonant and off-resonant pulses produce a
plasma core of roughly equal radius. The most important difference between the
two cases is the plasma sheath width $t_{0}$, which is much smaller for
resonant pulses (Fig. 7). Though this quantity fluctuates a lot more than
$r_{0}$, especially for the off-resonant case, it is clear that it is at least
2-3 times as large for off-resonant pulses.

Figure 6: Plasma core radius $r_{0}$ for a) $\mathcal{N}=5\cdot
10^{14}\mathrm{~{}cm^{-3}}$ vapor density and b) $\mathcal{N}=7\cdot
10^{14}\mathrm{~{}cm^{-3}}$ vapor density (left vertical axis). Solid red line
/ symbols mark resonant, blue dashed line / symbols mark off-resonant
measurements. Transmitted pulse D4$\sigma$ width is also shown for reference
with dotted lines (right vertical axis, red / blue for resonant / off-
resonant).  Figure 7: Plasma sheath width $t_{0}$ for a) $\mathcal{N}=5\cdot
10^{14}\mathrm{~{}cm^{-3}}$ vapor density and b) $\mathcal{N}=7\cdot
10^{14}\mathrm{~{}cm^{-3}}$. Red solid lines correspond to resonant, dashed
blue lines correspond to off-resonant pulses. Figure 8: Pulse energy ratio
$\mathcal{R}_{e}$ within the $r_{0}$ radius plasma column core for a)
$\mathcal{N}=5\cdot 10^{14}\mathrm{~{}cm^{-3}}$ vapor density and b)
$\mathcal{N}=7\cdot 10^{14}\mathrm{~{}cm^{-3}}$. Red solid lines correspond to
resonant, dashed blue lines correspond to off-resonant pulses.

Figure 6 shows that plasma is detected for some shots already at the start of
the breakthrough domain. These shots are rare initially, but their frequency
increases as $E_{in}$ progresses towards the confined beam domain. The
D4$\sigma$ width corresponding to the shots with plasma is small, whereas the
shots where plasma is not detected, are much wider. As $E_{in}$ increases,
shots with narrow D4$\sigma$ width and measurable plasma become more frequent
and the average transmitted pulse width drops quickly. By the end of the
breakthrough domain, all shots produce a measurable plasma column at the end
of the vapor source. Note that because our evaluation of the schlieren images
yields no direct information on the ionization fraction $P_{max}$, we only
assume that we have $P_{max}\approx 1$ at the plasma core if $r_{0}$ is
substantial (i.e. we have a plateau of the ionization fraction in the center).
Naturally, ’substantial’ must be in comparison with $t_{0}$, i.e we may safely
assume $P_{max}=1$, only if $r_{0}\gtrsim t_{0}$. This regime is reached much
sooner with resonant pulses (also) because of the smaller $t_{0}$. In fact for
resonant pulses we definitely reach it by the start of the confined beam
domain.

To estimate the significance of a thinner sheath, we may consider the
ionization probability profile $P_{ion}(r)$ described by Eqn. 1 integrated in
a plane perpendicular to the propagation direction. The energy loss of the
propagating pulse due to ionization per unit distance $dE^{ion}_{loss}/dz$
will be proportional to this quantity:

$\int
P_{ion}(r)rdrd\phi=r_{0}^{2}\pi+t_{0}^{2}\pi+\pi\sqrt{\pi}r_{0}t_{0}\sim\frac{dE^{ion}_{loss}}{dz}$
(7)

Here we assume that every ionization event removes exactly three photons from
the field of the laser pulse — the smallest number required for ionization at
these wavelengths. (The effects of above-threshold ionization with four or
more photons are thus excluded from this simple consideration, as is the
energy loss by other means e.g. atoms not ionized but left in an excited
electronic state.) Using values for $r_{0}$ and $t_{0}$ from Figs. 6 and 7
that correspond to $E_{in}=100$ mJ and inserting in Eq. 7, we can readily see
that $dE^{ion}_{loss}/dz$ is roughly twice as large for the off-resonant pulse
at this point. For lower (higher) pulse energies, where the ratio
$t_{0}/r_{0}$ is higher (lower), the ratio by which $dE^{ion}_{loss}/dz$ is
higher for off-resonant pulses will be higher (lower).

An interesting question regarding pulse propagation is the actual ratio of the
pulse energy that propagates inside the plasma column core. In traditional
filamentation in atmospheric gases, most of the energy propagates in the low
intensity wings of the propagating pulses around the plasma as the medium is
transparent for low-intensity light. In the resonant pulse scenario discussed
here, however, there is absorption for an arbitrarily low intensity and theory
predicts that most of the energy is channeled inside the core, where resonant
absorption is stopped by the removal of the valence electron. The actual ratio
$\mathcal{R}_{E}$ of pulse energy propagating within the plasma column core
can be estimated from the camera images by integrating the fluence
distribution within an $r_{0}$ radius around the center. Figure 8 shows the
results of this evaluation. For the resonant case, around 90 % of the pulse
energy is found to be channeled in the core in the asymptotic transparency
domain. Off-resonant pulses are channeled less efficiently, especially for low
energies. Nevertheless, even they show a behavior similar to resonant pulses
at high energy, i.e. they resemble the propagation of resonant pulses much
more than they do traditional non-resonant filamentation phenomena.

The experimental data thus proves, that resonant self focusing helps contain
laser pulse energy near the plasma column core very effectively. Pulse energy
is channeled in the core where the optical nonlinearity and absorption are
saturated due to complete single electron ionization. As a consequence, the
plasma sheath layer with partial ionization is much thinner for resonant
pulses. Because less energy is lost by resonant pulses for the partial
ionization, these pulses have a greater penetration depth for a given energy,
i.e. they require less energy to create a continuous plasma column of given
length. This is clearly very advantageous for creating a long plasma column
and particularly so for accelerator applications.

###  3.2 Considerations for longer plasma columns

Using the data we may also estimate the pulse energy that could be capable of
creating a $z_{max}=20$ m long plasma column. The estimate is based on the
observation made in numerical modeling of the system that at a given vapor
density and beam focusing, the properties of the propagating pulse (and thus
the properties of the plasma column) at a given spatial position will depend
only on the energy still remaining in the pulse at that point [15, 16]. It is
made by recursion using the information on Figs. 5 and 6 – the curves required
for the process in this example are shown together in Fig. 9 for convenience.
If we require a $r_{0}=0.5$ mm plasma column after $z_{max}=20$ m in
$\mathcal{N}=5\cdot 10^{14}\mathrm{~{}cm}^{-3}$ density vapor and consider
resonant ionizing pulses, we first establish that we need $E_{in}=29$ mJ to
have $r_{0}=0.5$ mm after $z=10$ m (left panel of Fig. 9). We now estimate
that if we have a pulse which still contains 29 mJ energy after $z=10$ m, that
pulse will be able to create the required plasma column in the second 10 meter
section of the 20-meter-long vapor as well. Finding $E_{in}^{\prime}$ that
yields $E_{out}=29$ mJ from the $E_{out}(E_{in})$ curve we obtain
$E_{in}^{\prime}=90$ mJ (right panel of Fig. 9, arrows with dotted line depict
the process). Thus, we estimate that a resonant pulse with $E_{in}\gtrapprox
90$ mJ would be able to create the required plasma column in a 20 m long vapor
source — a value well within the capabilities of the laser system.

Figure 9: Illustration of the recursive method for estimating the energy
necessary to create a 20 m long plasma column of a given radius in
$\mathcal{N}=5\cdot 10^{14}\mathrm{~{}cm}^{-3}$ density vapor. Left panel:
$r_{0}(E_{in})$ curve for resonant pulses from Fig. 6 a), rotated
anticlockwise by 90 degrees. Right panel: $E_{out}(E_{in})$ curve for resonant
pulses from Fig. 5 c). Dashed line in the right panel marks the line of best
fit for the upper range of the experimental data.

The same process can be used to estimate the energy required for a
$z_{max}=20$ m long, $r_{0}=1$ mm plasma column at the same vapor density.
Figure 9 left panel shows that $E_{in}=70$ mJ is required for $r_{0}=1$ mm
after $z=10$ m. Projecting this value as $E_{out}$ onto the right panel we see
that we are somewhat above the range of available experimental data. However,
using linear extrapolation, we can estimate that $E_{in}^{\prime\prime}\approx
152$ mJ yields $E_{out}=70$ mJ after 10 meters of propagation, which in turn
is the energy requirement for a pulse to create a 20 m long, 1 mm radius
plasma column.

Figure 10: Illustration of the recursive method for estimating the energy
necessary to create a 20 m long plasma column of a given radius in
$\mathcal{N}=5\cdot 10^{14}\mathrm{~{}cm}^{-3}$ density vapor for off-resonant
pulses. Left panel: $r_{0}(E_{in})$ curve for off-resonant pulses from Fig. 6
b), rotated anticlockwise by 90 degrees. Right panel: $E_{out}(E_{in})$ curve
for off-resonant pulses from Fig. 5 c). Dashed line in the right panel marks
the line of best fit for the upper range of the experimental data.

Repeating the process to investigate the possibility of generating a
$z_{max}=20$ m, $r_{0}\geq 0.5$ mm plasma column by off-resonant pulses, we
find that $E_{out}=56$ mJ is required after $z=10$ m propagation (Fig. 10).
However, projecting this value to find $E_{in}^{\prime}$ we are far above the
experimental curve due to the fact that the plasma column appears for larger
$E_{in}$ and transmitted energies are smaller for the off-resonant case.
Linear extrapolation again yields an estimate for the required pulse energy to
be $E_{in}^{\prime}\approx 157$ mJ, much greater than the 90 mJ value of the
resonant case. Also, this estimate is to be handled with greater caution due
to its distance from the measured range. The same process with the data for
$\mathcal{N}=7\cdot 10^{14}\mathrm{~{}cm}^{-3}$ density vapor and resonant
pulses readily yields $E_{in}^{\prime}\approx 118$ mJ for the limit of
creating a $z_{max}=20$ m plasma column with $r_{0}\geq 0.5$ mm (without any
extrapolation). We also note that given enough experimental data, the
recursive process could be repeated to estimate the pulse energy requirement
for 30 meter long or even longer plasma columns.

##  4 Simulation

###  4.1 Theoretical framework

A theory has been developed recently to describe the propagation of ultra-
short, ionizing laser pulses through Rb vapor under the condition of single
photon resonance from the ground state [15]. Results from numerical
simulations have been found to agree qualitatively with experimental findings
of resonant pulse propagation for the transmitted pulse energy and width [16].
Using this as a starting point, we have derived a more general theory to treat
the propagation of off-resonant pulses and resonant pulses in a unified way
and explore the difference between their behavior. Here we present only a very
brief outline, as the basic concept is the same as the one described in [16]
in greater detail.

The basic equation used for the pulse propagating along the $z$ direction is
written for the complex envelope function $\mathcal{E}(r,z,t)$ of the
axisymmetric laser field
$E(r,z,t)=\frac{1}{2}\mathcal{E}(r,z,t)\exp(ik_{0}z-\omega_{0}t)+c.c.$
($\mathcal{E}(r,z,t)$ ($\omega_{0}$ and $k_{0}$ being the central frequency
and wavenumber in vacuum). Standard methods in the treatment of ultrashort
pulses are employed [39] such as the paraxial approximation, transforming to
Fourier space with respect to $\tau=t-z/c$, the delayed time:
$\tilde{\mathcal{E}}(r,z,\omega)=\mathfrak{F}\\{\mathcal{E}(r,z,\tau)\\}$ and
using the Slowly Evolving Wave Approximation (SEWA) [40] to arrive at the
propagation equation:

$\displaystyle\partial_{z}\tilde{\mathcal{E}}=$
$\displaystyle\frac{i}{2k}\nabla_{\perp}^{2}\tilde{\mathcal{E}}+i\frac{k}{2\epsilon_{0}}\tilde{\mathcal{P}}$
(8)
$\displaystyle-\eta_{0}\hbar\omega_{0}\mathcal{N}\tilde{\mathcal{Q}}-\frac{ik}{2}\frac{e^{2}\mathcal{N}}{\epsilon_{0}m_{e}(\omega_{0}+\omega)^{2}}\tilde{\mathcal{R}}$

Here $e,m_{e}$ are the elementary charge and electron mass, $\epsilon_{0}$ is
the vacuum permittivity and $\eta_{0}$ the vacuum impedance. The first term on
the right-handside (RHS) of Eq. 8 is due to diffraction, the last term is the
plasma dispersion term. The third term on the RHS of Eq. 8 is an energy loss
term derived from the requirement that when an atom is ionized, an appropriate
number of photons are absorbed from the field (the number depending on the
electronic level that the valence electron was in prior to ionization).

Figure 11: a) Electronic levels of the rubidium atom that are included in the
theoretical model. The states used in the original, 4-level model [15] are
highlighted in red. Black arrows mark allowed dipole transitions between the
states, but transitions to different fine-structure sublevels are not
resolved, i.e. only a single arrow corresponds e.g. to the 780 nm D2 and 795
nm D1 lines from $\mathrm{5s^{2}S_{1/2}}$ to $\mathrm{5p^{2}P_{3/2}}$, and
$\mathrm{5p^{2}P_{1/2}}$ states.

The second term on the RHS is the atomic polarization, which, for the resonant
case, is dominated by Rabi-oscillation type transitions on single-photon
resonances. Contributions of this type cannot be expressed in terms of usual
optical nonlinear coefficients [41]. Therefore our theory includes an explicit
calculation of the atomic state using the Schrödinger equation written for the
$\alpha_{j}(t)$ probability amplitudes:
$|\psi\rangle=\sum_{j}\alpha_{j}(t)\exp(-i\omega_{j}t)|j\rangle$
($\hbar\omega_{j}$ is the energy of the energy eigenstate $|j\rangle$). The
original model [15] included just four atomic states, the ground state and
three excited states that are accessible via resonant interaction within the
bandwidth of the 780 nm ionizing pulses (see Figs. 2 and 11). However, the
dominance of these over other atomic states will be much less significant for
off-resonant pulses centered around 810 nm. Therefore the model has been
expanded to include 10 atomic levels in all (Fig. 11). The levels were chosen
by first considering a more general model with 18 atomic states and then
selecting only those which proved to acquire a maximum occupation probability
of at least 0.01 during the interaction with the laser. The evolution of the
atomic state is thus described by:

$\displaystyle\partial_{t}\alpha_{j}=$
$\displaystyle\frac{i}{2\hbar}\sum_{k}\mathcal{E}e^{-i\Delta_{jk}t}d_{jk}\alpha_{k}$
(9)
$\displaystyle+\frac{i}{2\hbar}\sum_{k^{\prime}}\mathcal{E}^{*}e^{i\Delta_{jk^{\prime}}t}d_{jk^{\prime}}\alpha_{k^{\prime}}-\frac{\Gamma_{j}}{2}\alpha_{j}$

The summation for index $k$ runs over the lower lying atomic states
($\hbar\omega_{k}<\hbar\omega_{j}$) for which the dipole matrix element
$d_{jk}\neq 0$, the summation for $k^{\prime}$ for higher lying states in a
similar manner. $\Delta_{jk}=\omega_{0}-(\omega_{j}-\omega_{k})$ is the
detuning of the laser central frequency from the
$|k\rangle\rightarrow|j\rangle$ transition frequency. Material parameters
needed in the calculation (level energies and dipole oscillator strengths)
have been obtained from the literature [42, 43, 44]. $\Gamma_{j}$ are level
loss rates due to ionization (intensity dependent), obtained from the so-
called PPT-formulas [45, 46, 47] and experimental data [48]. With the atomic
state evolution calculated, the atomic polarization term for Eq. 8 is given by
the dipole operator expectation value $\langle\hat{d}\rangle$:

$\tilde{\mathcal{P}}=\mathfrak{F}\left\\{\mathcal{N}\sum_{kl}\alpha_{k}^{*}\alpha_{l}d_{kl}\right\\}.$
(10)

where $\mathfrak{F}\\{\ldotp\\}$ marks the time Fourier transform.

###  4.2 Numerical results and discussion

We performed calculations for a pulse propagation scenario similar to the
experiments, with both resonant and off-resonant pulses. The computations were
done for a range of input pulse energies with vapor density
$\mathcal{N}=7\cdot 10^{14}\mathrm{~{}cm^{-3}}$. The input laser pulse was a
Gaussian beam, with focal parameters derived from the measured virtual laser
line camera fluence distributions as in [16]. A hyperbolic secant temporal
dependence of the initial field was assumed — note that this ideal pulse shape
leads to a significantly narrower spectral width for the pulse than that
actually measured (Fig. 2). Transmitted pulse properties and ionization
profiles were plotted after 10 meters of propagation for comparison with
experimental results.

Figure 12: Simulation results for $\mathcal{N}=7\cdot
10^{14}\mathrm{~{}cm^{-3}}$ vapor density a) transmitted pulse energy, b)
transmitted pulse $D4\sigma$ width, c) plasma core radius $r_{0}$ and d)
plasma sheath width $t_{0}$. Solid red lines mark 780 nm resonant results,
solid blue lines mark 810 nm off-resonant results. Experimental results are
also plotted with dashed lines (without error bars) for an easy comparison.
Experimental $t_{0}$ values on d) are rescaled by a factor of 1/3 to fit close
to the range of the simulation results.

Figure 12 a) and b) show transmitted pulse energies and widths, while Fig. 12
c) and d) shows plasma core radii and sheath widths. There is a good
qualitative similarity between experiment and theory albeit with some
quantitative discrepancy. Transmitted pulse energy (Fig. 12 a) ) is predicted
to be significantly larger for any given $E_{in}$ than that measured, but the
fact that the resonant pulse $E_{out}$ curve commences growing earlier and is
greater than its off-resonant counterpart is clearly represented in the
results. Transmitted pulse width (Fig. 12 b) ) also shows that simulation
reproduces well the breakthrough behavior of the pulses, the quick drop in
$D4\sigma$ to a minimum width, as well as the following expansion of the width
of the resonant pulse. The minimum width value obtained in the simulation is
somewhat smaller for the resonant case. The plasma $r_{0}$ behavior is also
well reproduced (Fig. 12 c) ), the curve for the resonant case displaying the
shoulder of near constant $r_{0}$ that corresponds approximately to a similar
region in the $D4\sigma$ width (the confined beam domain). Finally, the sheath
width $t_{0}$ (Fig. 12 d) ) that we obtained from simulation also shows that
the resonant pulse gives a $t_{0}$ value several times smaller than the off-
resonant pulse does. However, the numerical values from the simulation are
smaller by a factor of about three. (Experimental values have been rescaled on
Fig. 12 d), see also Fig. 7 b).)

Figure 13: Contour plot of the ionization probability at the vapor source exit
obtained from a), b) numerical simulation and c), d) schlieren imaging
measurements for $\mathcal{N}=7\cdot 10^{14}\mathrm{~{}cm^{-3}}$ density and
both resonant and off-resonant pulses.

It is noteworthy that the drop in transmitted pulse $D4\sigma$ is much sharper
for the simulation than measured experimentally. It has been shown in [16],
that the parameters of the ionizing pulses (measured beam waist parameters
$w_{0}$, $z_{0}$ and the beam profile shape as well) fluctuate somewhat shot-
to-shot. The simulations on the other hand are performed with constant beam
parameters. Varying the beam parameters in the simulation (within the range of
variation observed in the experiment) the sharp drop is smeared out and small-
scale features (present predominantly for the resonant pulse case) are
smoothed, averaged out [16].

To visualize the plasma column at the vapor source exit, we show a 2D plot of
the ionization probability (i.e. the plasma density) as a function of $E_{in}$
and the radius on Fig. 13. Plotted are: a), b) $P_{ion}(E_{in},r)$ at $z=10$ m
obtained from the simulation and c), d) ionization probabilities of the form
Eq. 1, with parameters $r_{0},t_{0}$ derived from bin averages of binned
$r_{0}$ and $t_{0}$ values obtained from the experiment (Figs. 6 and 7). The
plasma column core in the figures is the region close to the $r=0$ axis where
$P_{ion}=1$. The plots show good resemblance between simulation and experiment
— the fact that the experimental sheath layer width is much wider than the
simulated one is also visible. Finally, Figure 14 shows two plots for
simulated $P_{ion}(z,r)$ calculated for the entire length of the vapor source
in case of a resonant and an off-resonant pulse with the same energy. The
pulse energy is low, a little below the breakthrough transition of the
resonant pulse. The plots reflect very well that conclusions drawn in the
previous sections from quantities observed / calculated for the vapor source
exit hold during the entire propagation except for a small transient region
just after the entry. They depict explicitly that the resonant pulse of equal
energy creates a longer plasma column than the off-resonant one.

Figure 14: Ionization probability inside the vapor source from simulation
calculated for a $E_{in}=18$ mJ energy, a) 780 nm resonant and b) 810 nm off-
resonant pulse.

Finally, some comments on the possible causes of the large quantitative
discrepancy between numerically predicted and experimentally observed
$E_{out}$ and $t_{0}$ values. In several respects, the simulation assumes a
situation that is only imperfectly realized in the experiment. First, the real
laser pulse is not axially symmetric and its spatial distribution is far from
being a pure Gaussian beam. While non-axisymmetric (i.e. full 3D) simulations
are far out of scope, similar calculations have been conducted by us
previously with flattened Gaussian beams of various orders [49] using the
previous numerical model for resonant propagation [15].These studies
determined that the different transverse spatial dependence did not affect the
propagation substantially apart from a fairly short transition region at the
beginning of the vapor and did not result in a substantial difference in the
transmitted energy. Second, the idealized sech temporal pulse shape does lead
to a much narrower spectral width of the pulse than that measured for the
experimental input (Fig. 2). This may possibly cause a more significant
difference. Furthermore, the simplified calculation of ionization (required to
get an atomic model calculable in 2D propagation simulations) may also cause a
significant quantitative difference from the experiment.

##  5 Summary

We have presented both experimental and simulation results for the propagation
of resonant and off-resonant, ultra-short, ionizing laser pulses in rubidium
atomic vapor. The 780 nm central wavelength pulses were resonant with atomic
lines of the rubidium atom, while the 810 nm central wavelength off-resonant
pulses were just above the resonance wavelengths. We performed measurements
varying the input pulse energy from 0 to $\sim$150 mJ and repeated the
measurements for several vapor densities.

We measured the transmitted pulse energy and width after propagation along the
10-meter-long rubidium vapor column. Simultaneously, we have investigated the
transverse extent of the plasma column created by the pulses close to the
downstream end of the vapor column using schlieren imaging. From the schlieren
images, we determined the radius of the plasma column core, where the laser
pulses achieve one-electron ionization of the rubidium atoms with a
probability very close to 1. We also determined the plasma column sheath layer
width, which characterizes the fall-off distance of the plasma density from
the core to the unionized, neutral vapor.

We established that resonant pulses lose less energy during propagation per
unit distance and are confined more strongly in the transverse plane by the
interaction with the vapor. For low energy, resonant ionizing pulses create a
plasma with significantly larger core radius at the vapor’s end than off-
resonant pulses do. For high energy pulses, the radius of the plasma core is
approximately equal for the two cases. For the entire energy range
investigated, resonant pulses create a plasma column with a much thinner
sheath layer. Based on our results, we conclude that the greater energy loss
in case of off-resonant pulses comes predominantly from this wider sheath
layer, as they leave a much larger volume of partially ionized atoms around
the core. For this reason, resonant pulses are able to create a longer plasma
column than off-resonant pulses with the same energy. We observe that,
contrary to traditional filamentation in atmospheric gases, both resonant and
off-resonant pulses are channeled in the plasma core where further ionization
and thus energy absorption does not take place. However, resonant pulses are
being channeled much more efficiently with around 90% of the pulse energy
traveling in the core for mid- to high-energy pulses.

We have also used our experimental results to estimate the pulse energy
required for creating a 20-meter-long plasma column under identical
conditions. While the laser used in our experiments would be able to create
this double length plasma column with resonant pulses, it is far from it when
off-resonant pulses are used.

We generalized a theoretical description developed for the resonant pulse
propagation problem to treat the off-resonant case as well. Comparing
simulation results with measurement, we conclude that they are qualitatively
similar, even if the quantitative predictive power of the theory is still
lacking in certain respects. Using simulations results, we have shown that the
plasma column created by resonant pulses is bounded more sharply almost all
along the entire vapor column and that this allows the resonant pulses to
generate a much longer plasma column with equal initial pulse energy. Our
results show that single-photon resonances have a major effect on pulse
propagation and that they can be very advantageous if the requirement is to
generate a long, 10-meter-scale plasma column. Our results can potentially be
very significant for the construction of plasma wakefield accelerator devices,
but may also be of interest for other applications with high-power laser
pulses, e.g. lightning protection or remote sensing applications.

## Funding

The research was supported by the Hungarian National Research, Development and
Innovation Office (NKFIH) under the contract numbers NKFIH
2019-2.1.6-NEMZ_KI-2019-00004 and MEC_R 140947. On behalf of Project
Awakelaser we are grateful for the usage of ELKH Cloud [50] which helped us
achieve the results published in this paper.

## CRediT authorship contribution statements

G. Demeter: Conceptualization; Investigation; Formal analysis; Methodology;
Software; Data Curation; Visualization; Writing - original draft; J. T. Moody:
Conceptualization; Methodology; Resources; Supervision; Investigation; Data
Curation; Writing - review & editing; M. Á. Kedves: Conceptualization;
Investigation; Data Curation; Writing - review & editing; F. Batsch:
Resources; Writing - review & editing; M. Bergamaschi: Resources; Writing -
review & editing; V. Fedosseev: Resources; Writing - review & editing; E.
Granados: Resources; Writing - review & editing; P. Muggli: Project
Administration; Funding acquisition; Writing - review & editing; H. Panuganti:
Resources; Writing - review & editing; G. Zevi Della Porta: Resources; Writing
- review & editing;

## Declaration of competing interests

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * Couairon and Mysyrowicz [2007] A. Couairon and A. Mysyrowicz, Physics Reports 441, 47– (2007).
  * Bergé  _et al._ [2007] L. Bergé, S. Skupin, R. Nuter, J. Kasparian, and J.-P. Wolf, Reports on Progress in Physics 70, 1633– (2007).
  * Kandidov  _et al._ [2009] V. P. Kandidov, S A Shlenov, and O. G. Kosareva, Quantum Electronics 39, 205– (2009).
  * Joshi and Katsouleas [2003] C. Joshi and T. Katsouleas, Physics Today 56, 47 (2003).
  * Leemans and Esarey [2009] W. Leemans and E. Esarey, Physics Today 62, 44 (2009).
  * Malka  _et al._ [2008] V. Malka, J. Faure, Y. A. Gauduel, E. Lefebvre, A. Rousse, and K. T. Phuoc, Nature Physics 4, 447 (2008).
  * Esarey  _et al._ [2009] E. Esarey, C. B. Schroeder, and W. P. Leemans, Rev. Mod. Phys. 81, 1229 (2009).
  * Tajima  _et al._ [2020] T. Tajima, X. Q. Yan, and T. Ebisuzaki, Reviews of Modern Plasma Physics 4, 7 (2020).
  * Albert and Thomas [2016] F. Albert and A. G. R. Thomas, Plasma Physics and Controlled Fusion 58, 103001 (2016).
  * Hidding  _et al._ [2019] B. Hidding, B. Foster, M. J. Hogan, P. Muggli, and J. B. Rosenzweig, Phil. Trans. R. Soc. A. 377, 20190215 (2019).
  * Adli [2019] E. Adli, Phil. Trans. R. Soc. A. 377, 20180419 (2019).
  * Adli [2022] E. Adli, Journal of Instrumentation 17, T05006 (2022).
  * Gschwendtner  _et al._ [2016] E. Gschwendtner _et al._ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 829, 76– (2016), 2nd European Advanced Accelerator Concepts Workshop - EAAC 2015.
  * Adli  _et al._ [2018] E. Adli _et al._ (AWAKE Collaboration), Nature 561, 363 (2018).
  * Demeter [2019] G. Demeter, Phys. Rev. A 99, 063423 (2019).
  * Demeter  _et al._ [2021] G. Demeter, J. T. Moody, M. A. Kedves, B. Ráczkevi, M. Aladi, A.-M. Bachmann, F. Batsch, F. Braunmüller, G. P. Djotyan, V. Fedosseev, F. Friebel, S. Gessner, E. Granados, E. Guran, M. Hüther, V. Lee, M. Martyanov, P. Muggli, E. Öz, H. Panuganti, L. Verra, and G. Z. Della Porta, Phys. Rev. A 104, 033506 (2021).
  * Miao  _et al._ [2020] B. Miao, L. Feder, J. E. Shrock, A. Goffin, and H. M. Milchberg, Phys. Rev. Lett. 125, 074801 (2020).
  * Gupta  _et al._ [2022] D. N. Gupta, S. R. Yoffe, A. Jain, B. Ersfeld, and D. A. Jaroszynski, Scientific Reports 12, 20368 (2022).
  * Houard  _et al._ [2023] A. Houard, P. Walsch, T. Produit, _et al._ , Nature Photonics (2023), https://doi.org/10.1038/s41566-022-01139-z.
  * Muggli  _et al._ [2017] P. Muggli, E. Adli, R. Apsimon, _et al._ , Plasma Physics and Controlled Fusion 60, 014046 (2017).
  * Gschwendtner  _et al._ [2022] E. Gschwendtner, K. Lotov, P. Muggli, M. Wing, _et al._ (AWAKE Collaboration), Symmetry 14 (2022), 10.3390/sym14081680.
  * Öz and Muggli [2014] E. Öz and P. Muggli, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 740, 197 (2014), proceedings of the first European Advanced Accelerator Concepts Workshop 2013.
  * Plyushchev  _et al._ [2017] G. Plyushchev, R. Kersevan, A. Petrenko, and P. Muggli, Journal of Physics D: Applied Physics 51, 025203 (2017).
  * Batsch  _et al._ [2018] F. Batsch, M. Martyanov, E. Oez, J. Moody, E. Gschwendtner, A. Caldwell, and P. Muggli, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 909, 359 (2018), 3rd European Advanced Accelerator Concepts workshop (EAAC2017).
  * Settles [2001] G. S. Settles,  _Schlieren and Shadowgraph Techniques - Visualizing Phenomena In Transparent Media_ (Springer-Verlag, Berlin, Heidelberg, 2001).
  * Iwase  _et al._ [1998] O. Iwase, W. Süß, D. H. H. Hoffmann, M. Roth, C. Stöckl, M. Geissel, W. Seelig, and R. Bock, Physica Scripta 58, 634 (1998).
  * Clayton  _et al._ [1998] C. E. Clayton, K.-C. Tzeng, D. Gordon, P. Muggli, W. B. Mori, C. Joshi, V. Malka, Z. Najmudin, A. Modena, D. Neely, and A. E. Dangor, Phys. Rev. Lett. 81, 100 (1998).
  * Honda  _et al._ [2000] H. Honda, K. Katsura, E. Takahashi, and K. Kondo, Applied Physics B 70, 395 (2000).
  * Veloso  _et al._ [2006] F. Veloso, H. Chuaqui, R. Aliaga-Rossel, M. Favre, I. H. Mitchell, and E. Wyndham, Review of Scientific Instruments 77, 063506 (2006).
  * Batani  _et al._ [2019] D. Batani, J. Santos, P. Forestier-Colleoni, D. Mancelli, M. Ehret, J. Trela, A. Morace, K. Jakubowska, L. Antonelli, D. del Sorbo, M. Manclossi, and M. Veltcheva, Journal of Fusion Energy 38, 299 (2019).
  * Bachmann  _et al._ [2018] A.-M. Bachmann, M. Martyanov, J. Moody, A. Pukhov, and P. Muggli, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 909, 387–390 (2018), 3rd European Advanced Accelerator Concepts workshop (EAAC2017).
  * Kerscher [2021] M. Kerscher,  _Determination of Radial Plasma Column Profile from Examination of Diffraction Effects on Schlieren Imaging at AWAKE_ , Master’s thesis, Munich School of Engineering, Technical University of Munich (2021).
  * Saleh and Teich [2019] B. E. A. Saleh and M. C. Teich,  _Fundamentals of Photonics_ (John Wiley and Sons, Inc., 2019).
  * Siddons  _et al._ [2008] P. Siddons, C. S. Adams, C. Ge, and I. G. Hughes, Journal of Physics B: Atomic, Molecular and Optical Physics 41, 155004 (2008).
  * van Lange  _et al._ [2020] A. J. van Lange, P. van der Straten, and D. van Oosten, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 125402 (2020).
  * Adli  _et al._ [2019] E. Adli _et al._ (AWAKE Collaboration), Phys. Rev. Lett. 122, 054802 (2019).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * Bíró  _et al._ [2023] G. Bíró, M. Pocsai, I. F. Barna, G. G. Barnaföldi, J. T. Moody, and G. Demeter, Optics and Laser Technology 159, 108948 (2023).
  * Couairon  _et al._ [2011] A. Couairon, E. Brambilla, T. Corti, D. Majus, O. de J. Ramírez-Góngora, and M. Kolesik, The European Physical Journal Special Topics 199, 5– (2011).
  * Brabec and Krausz [1997] T. Brabec and F. Krausz, Phys. Rev. Lett. 78, 3282– (1997).
  * Boyd [2003] R. Boyd,  _Nonlinear Optics_ (Elsevier Science, 2003).
  * Steck [2009] D. A. Steck,  _Rubidium 85 D Line Data_ (available online at http://steck.us/alkalidata, (revision 2.1.2, 12 August 2009)).
  * Safronova  _et al._ [2004] M. S. Safronova, C. J. Williams, and C. W. Clark, Phys. Rev. A 69, 022509 (2004).
  * Kramida  _et al._ [2018] A. Kramida, Y. Ralchenko, J. Reader, and N. A. Team,  _NIST Atomic Spectra Database_ (National Institute of Standards and Technology, Gaithersburg, MD., 2018).
  * Perelomov  _et al._ [1966] A. M. Perelomov, V. S. Popov, and M. V. Terent’ev, Soviet Physics JETP 23, 924 (1966).
  * Perelomov  _et al._ [1967] A. M. Perelomov, V. S. Popov, and M. V. Terent’ev, Soviet Physics JETP 24, 207 (1967).
  * Perelomov and Popov [1967] A. M. Perelomov and V. S. Popov, Soviet Physics JETP 25, 336 (1967).
  * Duncan  _et al._ [2001] B. C. Duncan, V. Sanchez-Villicana, P. L. Gould, and H. R. Sadeghpour, Phys. Rev. A 63, 043411 (2001).
  * Gori [1994] F. Gori, Optics Communications 107, 335 (1994).
  * Héder  _et al._ [2022] M. Héder, E. Rigó, D. Medgyesi, R. Lovas, S. Tenczer, F. Török, A. Farkas, M. Emődi, J. Kadlecsik, G. Mező, Á. Pintér, and P. Kacsuk, Információs Társadalom 22, 128 (2022).