Nonlinear effects in photoionization over a broad photon-energy range within the TDCIS scheme

TL;DR

This paper reviews the TDCIS method for modeling nonlinear photoionization across a broad photon-energy spectrum, highlighting its application to argon and xenon, and exploring phenomena like giant resonances and multiphoton ionization in the x-ray regime.

Contribution

It provides a comprehensive tutorial on the TDCIS scheme's application to nonlinear ionization, including wave function propagation, photoelectron spectra calculation, and analysis of high-energy phenomena.

Findings

Analysis of xenon 4d giant dipole resonance reveals resonance substructure.

Identification of the intensity threshold for multiphoton ionization at hard x-ray energies.

Application of TDCIS demonstrates its effectiveness across a broad photon-energy range.

Abstract

The present tutorial provides an overview of the time-dependent configuration interaction singles (\acs{TDCIS}) scheme applied to nonlinear ionization over a broad photon-energy range. The efficient propagation of the wave function and the calculation of photoelectron spectra within this approach are described and demonstrated in various applications. Above-threshold ionization of argon and xenon in the extreme ultraviolet energy range is investigated as an example. A particular focus is put on the xenon $4d$ giant dipole resonance and the information that nonlinear ionization can provide about resonance substructure. Furthermore, above-threshold ionization is studied in the x-ray regime and the intensity regime, at which multiphoton ionization starts to play a role at hard x-ray photon energies, is identified.

Tables

Table 1. Table 1: Xenon resonances in the intra- and interchannel cases, data taken from Ref.  Chen et al. ( 2015 ) .

scheme	resonance	Energy (eV)	“lifetime” (as)
intrachannel	$R_1,R_2,R_3$	77	60
interchannel	$R_1$	74.3	26
	$R_2$	107.6	11

Table 2. Table 2: Two-photon ATI cross sections for the K 𝐾 K shell of the light elements (E is the photon energy). Data reproduced from Ref.  Tilley et al. ( 2015 ) with permission. All were calculated from integrating the corresponding photoelectron peaks at 10 20 superscript 10 20 10^{20}  Wcm -2 intensity and 0.12 0.12 0.12  fs of pulse duration.

E (keV)	$K$-shell generalized two-photon cross sections (cm4s)
	Hydrogen	Carbon	Nitrogen	Oxygen
8	$1.44\times {10}^{-66}$	$1.64\times {10}^{-62}$	$3.21\times {10}^{-62}$	$5.62\times {10}^{-62}$
10	$4.69\times {10}^{-67}$	$4.61\times {10}^{-63}$	$9.23\times {10}^{-63}$	$1.70\times {10}^{-62}$
12	$1.72\times {10}^{-67}$	$1.79\times {10}^{-63}$	$3.82\times {10}^{-63}$	$6.94\times {10}^{-63}$

Equations

$$M_{F\leftarrow I}^{(2)}=\sum _K\frac{⟨F|{\hat{H}}_{\mathrm{int}}|K⟩⟨K|{\hat{H}}_{\mathrm{int}}|I⟩}{E-E_K+\frac{i}{2}{\Gamma }_K+E_I},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{P}_0(t)=-P_N(t)=-{\sigma }^{(N)}{j}^N{P}_0(t),$$

$$P_0(t)=\exp \left(-{\int }_{-\infty }^t\mathrm{d}\tau {\sigma }^{(N)}{j}^N\right)+const.$$

$$\ln P_N=N\ln I+\ln {\sigma }^{(N)}+const.$$

$$i\frac{\partial }{\partial t}|{\Psi }^N(t)⟩=\hat{H}(t)|{\Psi }^N(t)⟩,$$

$$\hat{H}(t)=\underset{{\hat{H}}_0}{\underbrace{\sum _{n=1}^N\left(\frac{{\hat{\mathbf{p}}}_n^2}{2}-\frac{Z}{|{\hat{\mathbf{r}}}_n|}+{\hat{V}}_{\mathrm{MF}}({\hat{\mathbf{r}}}_n)\right)}}+\underset{{\hat{H}}_1}{\underbrace{\frac{1}{2}\sum _{n\ne n^{\prime }}\frac{1}{|{\hat{\mathbf{r}}}_n-{\hat{\mathbf{r}}}_{n^{\prime }}|}-\sum _{n=1}^N{\hat{V}}_{\mathrm{MF}}({\hat{\mathbf{r}}}_n)}}+\underset{{\hat{H}}_{\mathrm{int}}}{\underbrace{\hat{\mathbf{p}}\cdot \mathbf{A}(t)}},$$

$$|{\Phi }_0^N⟩=\det \left(|{\phi }_1⟩,\dots ,|{\phi }_N⟩\right),$$

$$|{\Phi }_i^a⟩={\hat{c}}_a^{†}{\hat{c}}_i|{\Phi }_0⟩.$$

$$|\Psi (t)⟩$$

$$i{\dot{\alpha }}_0(t)$$

$$i{\dot{\alpha }}_i^a(t)$$

$$+\mathbf{A}(t)\cdot \left(⟨{\Phi }_i^a|\hat{\mathbf{p}}|{\Phi }_0⟩{\alpha }_0(t)+\sum _{j,b}⟨{\Phi }_i^a|\hat{\mathbf{p}}|{\Phi }_j^b⟩{\alpha }_j^b(t)\right).$$

$$i{\dot{\alpha }}_0$$

$$i{\dot{\alpha }}_i^a$$

$$-A(t)\sum _{i^{\prime }}{p}_{i^{\prime }i}{\alpha }_{i^{\prime }}^a,$$

$$v_{pqrs}=\int {\mathrm{d}}^{3}r{\mathrm{d}}^3{r}^{\prime }{\phi }_p^{†}(\mathbf{r}){\phi }_q^{†}({\mathbf{r}}^{\prime })\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\phi }_r(\mathbf{r}){\phi }_s({\mathbf{r}}^{\prime }),$$

$$p_{ab}=⟨{\phi }_a|{\hat{p}}_z|{\phi }_b⟩.$$

$$V(\mathbf{r})=-\frac{Z}{r}+\int {\mathrm{d}}^3{r}^{\prime }\frac{\rho (\mathbf{r})}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+V_{\mathrm{ex}}(\mathbf{r}),$$

$$V(\mathbf{r})=-(Z^{\prime }+1)/r,\text{ if }-Z/r+\int {\mathrm{d}}^3{r}^{\prime }\frac{\rho (\mathbf{r})}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+V_{\mathrm{ex}}(\mathbf{r})<-(Z^{\prime }+1)/r,$$

$$R(\mathbf{x})=\frac{{\mathbf{x}}^{\ast }A\mathbf{x}}{{\mathbf{x}}^{\ast }\mathbf{x}},$$

$${\lambda }_1$$

$${\lambda }_n$$

$$\nabla R(\mathbf{x})=\frac{2}{{\mathbf{x}}^{\ast }\mathbf{x}}\left[A\mathbf{x}-R(\mathbf{x})\mathbf{x}\right],$$

$$\mathrm{span}\left\{{\mathbf{q}}_1,\dots ,{\mathbf{q}}_{k+1}\right\}=\mathrm{span}\left\{{\mathbf{q}}_1,A{\mathbf{q}}_1,A^2{\mathbf{q}}_1,\dots ,A^k{\mathbf{q}}_1\right\},$$

$$|\Psi (t)⟩=\hat{U}(t,0)|\Psi (0)⟩,$$

$$\hat{U}(t,0)=\mathcal{T}\exp \left[-i{\int }_0^t\mathrm{d}\tau \hat{H}(\tau )\right]$$

$$\mathcal{T}[\hat{H}(t_1)\hat{H}(t_2)]=\{\begin{array}{ll}\hat{H}(t_1)\hat{H}(t_2)& \text{mbox}if{t}_1<t_2\\ \hat{H}(t_2)\hat{H}(t_1)& \text{mbox}if{t}_1>t_2.\end{array}$$

$$|\Psi (t+\mathrm{d}t)⟩=e^{-i\hat{H}\left(t+\frac{\mathrm{d}t}{2}\right)\mathrm{d}t}|\Psi (t)⟩+\mathcal{O}(\mathrm{d}{t}^3).$$

$${\hat{H}}^n|\Psi (0)⟩={\mathbf{v}}^n,$$

$$Q^T\hat{H}Q$$

$$U{e}^{-iD\mathrm{d}t}{U}^T$$

$$e^{-i\hat{H}\mathrm{d}t}|\Psi ⟩$$

$$T=\left(\leavevmode\hbox to124.16pt{\vbox to115.63pt{\pgfpicture\makeatletter\hbox{\hskip 62.07799pt\lower-57.81326pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } {}{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{\offinterlineskip{}{}{{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}}{{{}}}{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-58.74498pt}{-54.48026pt}\pgfsys@invoke{ }\hbox{\vbox{\halign{\pgf@matrix@init@row\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding&&\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding\cr\hfil\hskip 7.93152pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.59851pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\alpha_{0}}$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip 7.93152pt\hfil&\hfil\hskip 10.06113pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-6.72812pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\ \beta_{1}}$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip 10.06113pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip 7.56113pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.22812pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\beta }_1$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 7.56113pt\hfil&\hfil\hskip 10.43152pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.09851pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\alpha }_1$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.43152pt\hfil&\hfil\hskip 15.06113pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-11.72812pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\beta }_2$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 15.06113pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 10.78304pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.45003pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${{\text{\Huge 0}}}$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.78304pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 10.06113pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-6.72812pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\beta }_2$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.06113pt\hfil&\hfil\hskip 15.43152pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-12.09851pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\alpha }_2$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 15.43152pt\hfil&\hfil\hskip 7.56113pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.22812pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\beta }_3$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 7.56113pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 11.04945pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.71645pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\beta }_{N-1}$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 11.04945pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip 10.78304pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.45003pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$\text{0}$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.78304pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 3.33301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 3.33301pt\hfil&\hfil\hskip 11.04945pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.71645pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\beta }_{N-1}$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 11.04945pt\hfil&\hfil\hskip 9.0865pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.7535pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${\alpha }_N$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 9.0865pt\hfil\cr}}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}{{{{}}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}{}{}}{}{ {}{}{}} {{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{ {}{}}{}{}{{}{}}}}}{{}}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setdash{0.4pt,4.0pt}{0.0pt}\pgfsys@invoke{ }{}\pgfsys@moveto{5.50359pt}{-9.78189pt}\pgfsys@lineto{41.2552pt}{-46.64171pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{ {}{}{}}{}{ {}{}{}} {{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{ {}{}}{}{}{{}{}}}}}{{}}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setdash{0.4pt,4.0pt}{0.0pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-18.40222pt}{-9.92189pt}\pgfsys@lineto{17.2562pt}{-44.0028pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{ {}{}{}}{}{ {}{}{}} {{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{ {}{}}{}{}{{}{}}}}}{{}}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setdash{0.4pt,4.0pt}{0.0pt}\pgfsys@invoke{ }{}\pgfsys@moveto{33.35774pt}{-9.74748pt}\pgfsys@lineto{39.76276pt}{-16.18791pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\right).$$

Full text

Acronyms

ADK Ammosov, Delone, and Krainov APS American Physical Society ARPACK Arnoldi package ATI Above-threshold ionization CAP Complex absorbing potential CI Configuration interaction CIS Configuration interaction singles ECS Exterior complex scaling FEL Free-electron laser FERMI Free-electron Laser for Multidisciplinary Investigations, Trieste FLASH Free-electron Laser Hamburg FWHM Full width at half maximum GDR Giant dipole resonance HFS Hartree-Fock-Slater HHG High harmonic generation IOP Institute of Physics LAPACK Linear algebra package LCLS (Stanford) LINAC Coherent Light Source MBES Magnetic bottle electron spectrometer OCT Optimal control theory PAD Photoelectron angular distribution PES Photoelectron spectrum RPAE Random-phase approximation with exchange SACLA Spring-8 Ångström Compact Free Electron Laser SAE Single-active electron SASE Self-amplified spontaneous emission SES Smooth exterior scaling SLAC Stanford Linear Accelerator Center SPring-8 Super Photon Ring 8 GeV t-surff Time-dependent surface flux TDCIS Time-dependent configuration interaction singles UV Ultraviolet VUV Vacuum Ultraviolet XATOM Integrated toolkit for x-ray atomic physics XCID Configuration interaction dynamics package for multichannel dynamics XFEL X-ray free-electron laser XUV Extreme ultraviolet

I Introduction

Remarkable progress in the realm of light-source development has shaped the physics of the 20th and the beginning of the 21st century. The power of producing coherent electromagnetic radiation of high intensity evolved simultaneously with the quest for ever increasing precision and for more information about the structure and dynamics of matter. Today, modern light sources deliver intense, ultrashort pulses at frequencies ranging from the terahertz to the hard x-ray regime. They provide the experimental means to control and image atomic and molecular systems and to test theoretical predictions of nonlinear processes Corkum and Krausz (2007); Chapman et al. (2006); Itatani et al. (2004); Spanner et al. (2004). On the one hand light is utilized to investigate the structure and dynamics of atoms; on the other hand it is employed to control atomic degrees of freedom and to steer electrons, which in turn leads to the development of new technology, such as attosecond light sources Giles (2002); Corkum and Krausz (2007); Krausz and Ivanov (2009). In the interaction with light the forces exerted on the electrons can be comparable to the intra-atomic forces and the ultrashort pulse durations reach the typical time scales involved in electronic excitations in atoms, molecules and clusters, roughly between 50 attoseconds (1\leavevmode\nobreak\as$\leavevmode\nobreak\ =10^{-18}$ s) and 50 femtoseconds (1\leavevmode\nobreak\fs $=10^{-15}$ s) Krausz and Ivanov (2009); Wabnitz et al. (2002). Therefore, the realm of strong-field physics, multiphoton processes and light-matter interactions on an ultrashort timescale, ranging from fs to as, has become a focus of interest. All of these processes are directly linked to the process of photo-excitation and -ionization.

Belonging to the last generation of light sources, free-electron lasers (FELs) provide extraordinarily intense light pulses that permit the investigation and control of inner-shell processes, Auger decay or above-threshold ionization (ATI) Kanter et al. (2011); Auger (1975); Agostini et al. (1979). Typically, the photon energies range from the UV to the x-ray range and the duration of FEL pulses are as short as a few femtoseconds, i.e., suitable for studying phenomena in atoms and molecules with a new quality in time resolution.

The present tutorial deals with the theoretical investigation of the nonlinear response of atomic systems interacting with intense light pulses, spanning a broad frequency range from the infrared to the x-ray regime. After an introduction to the theoretical and computational framework an analysis of the adiabatic eigenstates of the many-body Hamiltonian will be performed in order to clarify the notion of adiabaticity in strong-field ionization. Subsequently, the calculation of photoelectron spectra will be presented and applied to XUV ionization of argon. Xenon will also be investigated in the XUV energy regime, concentrating on the $4d$ giant dipole resonance (GDR). Finally, the x-ray regime will be addressed in order to quantify the impact of multiphoton ionization in current and future FEL experiments.

For the above-mentioned purposes, the theoretical framework must be capable of efficiently describing the absorption, emission and the scattering of photons by atoms. Consequently, a good description of the electronic structure, of the (strong) light field and of the interaction of light with matter is needed. Moreover, in the present article aspects of many-body physics in the atomic shell shall be investigated. Phenomena which involve collective electronic behavior cannot be described within the single-active electron picture. Correlation effects in the atomic shell can be scrutinized by employing the configuration- interaction singles (CIS) scheme. The wave function is expanded in the CIS basis and the time-dependent Schrödinger equation is solved “ab initio” by numerical time propagation. Cases where collective electronic behavior plays a crucial role will be exemplified for xenon.

Throughout this tutorial atomic units will be used, which are the “natural” units for processes occurring on atomic time and length scales. In atomic units the mass and the charge of the electron are set to unity, i.e., $m_{e}=1$, and $|e|=1$. Furthermore, the reduced Planck constant is set to one, $\hbar=1$. It follows that the fine-structure constant $\alpha=\frac{e^{2}}{\hbar c}=\frac{1}{c}\approx\frac{1}{137}$, where $c$ is the speed of light. The Bohr radius, or bohr, is the atomic unit of length, $a_{0}=\frac{1}{\alpha}\frac{\hbar}{m_{e}c}=0.53$ Å and the Hartree energy, or hartree, is the atomic unit of energy, $E_{h}=\frac{\hbar^{2}}{m_{e}a_{0}^{2}}=27.21$ eV. The atomic unit of time corresponds to the time the electron needs to travel the distance of a Bohr radius at the velocity $v_{0}=\frac{a_{0}E_{h}}{\hbar}$, such that $t_{0}=\frac{\hbar}{E_{h}}=\frac{1}{\alpha^{2}}\frac{\hbar}{m_{e}c^{2}}=24$ as. Then, the classical orbital period of the electron in the first Bohr orbit is given by $T=2\pi t_{0}\approx 150$ as.

II Photoionization

As one of the most probable processes to happen when light interacts with matter, photoionization has been studied extensively from the very beginning of quantum theory Bethe and Salpeter (1977). Both, theoretically and experimentally, the ionization of atoms and molecules has served as a tool to investigate the nature of the atomic shell and the molecular orbitals Starace (1982). If the light fields interacting with matter are strong, e.g. at FELs, multiphoton processes play a significant role  Young et al. (2010); Moshammer et al. (2007). The removal of a deep inner-shell electron is followed by various processes depending on the atomic states and the photon energy Drescher et al. (2002); Fukuzawa et al. (2013).

There are two ionization regimes which can be classified according to the so-called Keldysh parameter $\gamma=\sqrt{I_{p}/(2U_{p})}$ Ivanov et al. (2005): The “tunneling ionization” Keldysh (1965) and the “(perturbative) multiphoton ionization” Mainfray and Manus (1991). $I_{p}$ is the ionization potential and $U_{p}=I/(4\omega^{2})$ is the ponderomotive potential, which corresponds to the average energy of a free electron oscillating in the electric field. It involves the intensity of the light field, $I$, and the central frequency, $\omega$. In this way, the Keldysh parameter sets into proportion the frequency of the light field and the field amplitude, $\gamma=2\omega\sqrt{I_{p}/(2I)}$. For $\gamma\ll 1$ the ponderomotive potential is much larger than the ionization potential and ionization is governed by tunneling ionization, while for $\gamma\gg 1$ the process is governed by multiphoton ionization. In the range of $\gamma\approx 1$ both effects compete with each other Yamanouchi (2010).

Regarding the theoretical description there is a significant difference between the two regimes: In general, the tunneling regime requires a nonperturbative description, because the magnitude of the light’s electric field is comparable with the intra-atomic electric field and cannot be treated as a small perturbation. The electric field, which oscillates with a small frequency, bends the Coulomb potential in such a way that a barrier of finite width is created through which the electron can tunnel and leave the ion. This is shown schematically in Fig. 1 a). Since this is a quasistatic picture, the Keldysh parameter has been connected to the notion of adiabaticity in the ionization process Mevel et al. (1993). Far into the tunneling regime, the atomic response is considered to be purely adiabatic, which in this context means that the ionization rate at a given time is solely defined by the instantaneous electric field.

For relatively weak fields compared to the intra-atomic Coulomb potential (for atomic systems $I\ll 10^{16}$Wcm*-2*) the light field can be treated as a perturbation and the light-matter interaction can be classified within the lowest orders of perturbation theory. It can be viewed as a process where the system absorbs simultaneously several photons to undergo ionization, see Fig. 1 b). The order of the interaction or, equivalently, the number of photons that are absorbed is described by the corresponding order of perturbation theory. In this regime cross sections for the absorption or scattering of photons by atoms can be defined Saenz and Lambropoulos (1999). For this purpose, the transition matrix element $M_{F\leftarrow I}$ associated with the transition induced by the interaction Hamiltonian $\hat{H}_{\mathrm{int}}$ from an initial state of the coupled atom-light system $I$ to a final state $F$ is calculated. For instance, the absorption of two photons is described in second-order perturbation theory. All possible pathways leading to the final state by absorbing two photons have to be taken into account. Through the absorption of one photon an intermediate state $K$ is populated and a summation over all those possible intermediate states must be carried out

[TABLE]

where $E$ is the energy of the photon, $E_{K}$ and $\Gamma_{K}$ are the energy and the decay rate of the intermediate state, respectively, and $E_{I}$ is the energy of the initial state. $\Gamma_{K}$ is associated to the state’s life time $\tau_{K}$ by the relation $\Gamma_{K}=1/\tau_{K}$. The two-photon cross section, $\sigma^{(2)}$, is proportional to $|M^{(2)}_{F\leftarrow I}|^{2}$.

Exemplarily and for simplicity disregarding indirect channels the rate equation for the population evolving due to an $N$-photon ionization process reads

[TABLE]

which involves the ground state population $P_{0}$, the product of the cross section for the $N$-photon process, $\sigma^{(N)}$, and the $N^{\mathrm{th}}$ power of the photon flux, $j=I/\omega_{\mathrm{photon}}$, i.e., the number of incident photons per unit time and unit area. Usually, the $\sigma^{(N)}$ for $N>1$ are called “generalized” cross sections, since their units are not units of a cross section. The differential equation is solved through integration

[TABLE]

It follows that as long as the saturation regime is not reached, i.e., as long as the intensity is low enough that the ground state is not depleted (i.e., $P_{0}\approx 1$) and the cross section does not depend on the intensity, the following relation holds:

[TABLE]

The linear dependence results in a straight line on a double-logarithmic scale and allows to read off the order of the ionization process from the slope of the curve (i.e., one for one-photon processes, two for two-photon processes, and so forth).

III Theory and Method

The general time dependent Schrödinger equation of an $N$-electron system is given by

[TABLE]

where $|\Psi^{N}(t)\rangle$ is the $N$-electron wave function. Here, the Hamiltonian has the form

[TABLE]

where the single-electron part of the Hamiltonian, in the following denoted by $\hat{H}_{0}$, contains the kinetic energy $\hat{T}\equiv\sum_{n}\hat{\mathbf{p}}_{n}^{2}/2$, the nuclear potential $\hat{V}_{\mathrm{nuc}}\equiv-\sum_{n}Z/|\hat{\mathbf{r}}_{n}|$ and the potential at the mean-field level $\hat{V}_{\mathrm{MF}}$. The electron-electron Coulomb interaction $\hat{V}_{\mathrm{e-e}}\equiv\frac{1}{2}\sum_{n\neq n^{\prime}}1/|\hat{\mathbf{r}}_{n}-\hat{\mathbf{r}}_{n^{\prime}}|$ completes the many-electron part of the Hamiltonian $\hat{H}_{1}\equiv\hat{V}_{\mathrm{e-e}}-\hat{V}_{\mathrm{MF}}$. Finally, $\hat{H}_{\mathrm{int}}$ denotes the light-matter interaction in the minimal coupling scheme involving the momentum operator $\hat{\mathbf{p}}=\sum_{n}\hat{\mathbf{p}}_{n}$ and the vector potential ${\mathbf{A}}(t)$. The Coulomb gauge for the vector potential, $\mathrm{div}\mathbf{A}\!=\!0$, is often chosen in the context of atomic physics. As a consequence, the vector potential is divergence-free and purely transverse. Here, the light is described semiclassically and as a function of time only. Due to the exclusive time dependence the electric field does not change over the spatial extent of the atom and, consequently, our description remains in the dipole approximation.

III.1 Time-dependent CI singles (TDCIS) scheme

In order to solve the electronic structure problem a convenient basis set must be found in which the wave function can be expanded. In quantum chemistry a widely-used scheme is the configuration interaction (CI) Szabo and Ostlund (1996): Starting from the Hartree-Fock ground state $|\Phi^{\mathrm{HF}}_{0}\rangle$ the configuration space is built up by exciting electrons from the occupied orbital and promoting them to an orbital which was previously unoccupied (also called a virtual orbital). This yields a one-particle–one-hole state $|\Phi_{i}^{a}\rangle$ for the $N$-electron system. Exciting a second electron and letting it occupy another previously unoccupied orbital yields a two-particle–two-hole state $|\Phi_{i,j}^{a,b}\rangle$, promoting a third electron from an initially occupied orbital to a virtual orbital represents a three-particle–three-hole state $|\Phi_{i,j,k}^{a,b,c}\rangle$, and so forth. This is visualized schematically in Fig. 2, where $i,j,k,\ldots$ denote an initially occupied orbital and $a,b,c,\ldots$ denote a virtual orbital.

The configuration space is complete if all possible particle–hole excitations are taken into account. However, for higher excitations the space of virtual orbitals grows immensely. Therefore, neglecting higher order excitations $|\Phi_{i,j}^{a,b}\rangle,|\Phi_{i,j,k}^{a,b,c}\rangle,\ldots$, we truncate the configuration space and consider only states where one, single electron is excited to a previously unoccupied state. This space is called configuration interaction singles (CIS) space.

To construct it we start from the Hartree-Fock ground state $|\Phi_{0}^{N}\rangle$ of a closed-shell $N$-electron system as a reference state which is obtained from the vacuum state $|0\rangle$ by applying electron creation operators. The anticommuting operators $\hat{c}_{p\sigma}^{\dagger}$ and $\hat{c}_{p\sigma}$ are spin-orbital creation and annihilation operators, which create an electron in the spin orbital $|\varphi_{p\sigma}\rangle$ or annihilate an electron from this orbital, respectively, i.e., $\hat{c}^{\dagger}_{p\sigma}|0\rangle=|\varphi_{p\sigma}\rangle$. The field-free one-particle Hamiltonian $\hat{H}_{0}$ has the form $\hat{H}_{0}=\sum_{p}\varepsilon_{p}\sum_{\sigma}\hat{c}^{\dagger}_{p\sigma}\hat{c}_{p\sigma}$, such that $\hat{H}_{0}|\varphi_{p\sigma}\rangle=\varepsilon_{p}|\varphi_{p\sigma}\rangle$, where $\varepsilon_{p}$ denotes the energy of the orbital $|\varphi_{p\sigma}\rangle$. We consider closed-shell atoms and processes in which the total spin of the system is not altered ($S=0$, and there is no magnetic field involved), such that only spin singlets occur. Therefore, we henceforth drop the spin index and treat only the spatial part of the orbitals $|\varphi_{p}\rangle$. The Hartree-Fock ground state in our notation is the antisymmetrized product, or Slater determinant

[TABLE]

which involves the $N$ energetically lowest spin orbitals. In order to account for electronic excitations within the CIS scheme we build the so-called one-particle–one-hole configurations which have the form

[TABLE]

The index $i$ symbolizes an initially occupied orbital and $a$ denotes a virtual orbital in the sense described previously. The total N-electron wave function (now omitting the superscript $N$ for better legibility) for solving Eq. (5) is expanded within the CIS approach as Greenman et al. (2010)

[TABLE]

Inserting this expansion into the Schrödinger equation (5) and projecting it onto the states $|\Phi_{0}\rangle$ and $|\Phi_{i}^{a}\rangle$ yields the following equations of motion for the time-dependent expansion coefficients $\alpha_{i}^{a}(t)$ Greenman et al. (2010):

[TABLE]

The computational challenge lies in the numerical evaluation of the Coulomb matrix elements $\langle\Phi_{i}^{a}|\hat{H}_{1}|\Phi_{j}^{b}\rangle$, which must be calculated for all active occupied and virtual orbitals. In order to determine all the virtual orbitals that must be included in the calculation, a cut-off energy $E_{\mathrm{cut}}$ is defined up to which the virtual orbitals are calculated. Of course, it is essential to make the virtual space large enough to contain all states of interest.

In the following we assume that the electric field is linearly polarized along the $z$-axis, such that the light-atom interaction term simplifies to the projection of the momentum on the direction of the vector potential, $\hat{p}_{z}A_{z}$. Using the Slater-Condon rules Szabo and Ostlund (1996) and writing the one- and two-body matrix elements explicitly in terms of the spatial orbitals the Eqs. (10) read Rohringer et al. (2006):

[TABLE]

where the two-body matrix elements are given by

[TABLE]

and the matrix elements of the dipole operator, which is a one-body operator, are of the form

[TABLE]

The differential equations (11) are solved by numerical time propagation either using the $4^{\mathrm{th}}$ order Runge-Kutta algorithm or the Lanczos propagation, which will be introduced in the next section. Since the light-matter interaction is included in a nonperturbative manner in the Hamiltonian and, consequently, in Eq. (5) all orders of the interaction are automatically included in the calculations. In principle, this approach can cover all frequency ranges as long as the dipole approximation can be assumed to be valid.

III.2 Multichannel physics and electron correlations

Collective effects or correlation phenomena are abundant in atomic and molecular systems. Examples are autoionization, Auger decay, Fano resonances, the giant dipole resonance in xenon, and interatomic Coulombic decay Becker and Shirley (1996); Fano (1961); Fetter and Walecka (1980); Cederbaum et al. (1997). All of these involve the description of electronic dynamics beyond the single-active electron approach.

The CIS approach encapsulates many-body interactions through the coupling of different configurations contained in the Coulomb matrix elements $\langle\Phi_{i}^{a}|\hat{V}_{\mathrm{e-e}}|\Phi_{j}^{b}\rangle$. The set of active occupied orbitals $i$ builds the space of channels through which excitation and ionization can occur which automatically allows for the distinction and analysis of multichannel physics. Through the coefficients $\alpha_{i}^{a}(t)$ ionization and excitation processes can be analyzed in a channel-resolved manner. In this way, channel-resolved quantities can be inferred from the $N$-electron wave function, such as ionization probabilities, cross sections, and the ion density matrix (see, e.g., Ref. Pabst et al. (2011)).

In theory, we can manually switch off the correlation effects. If we allow only for the Coulomb matrix elements with hole indices $i=j$ to be nonzero, $\langle\Phi_{i}^{a}|\hat{V}_{\mathrm{e-e}}|\Phi_{i}^{b}\rangle$, we obtain the “intrachannel” picture, which is visualized in Fig. 3​ a). Physically, this means that the excited electron, although interacting with the hole, cannot modify the ionic hole state (depicted by open circles) in the remaining parent ion.

Without this restriction there can be non-vanishing matrix elements $\langle\Phi_{i}^{a}|\hat{V}_{\mathrm{e-e}}|\Phi_{j}^{b}\rangle$ for two configurations whose hole indices are different from one another, additionally to a different index of the excited electron, $i=j$ and $i\neq j$. This is called “interchannel” coupling and means that a simultaneous change of the excited electron state and the ionic hole state through electron-electron correlation is permitted. Fig. 3 ​b) depicts this situation, where states with $i\neq j$ and $a\neq b$, are coupled. In this way a correlated particle-hole pair is created. Whenever it is not possible to write the wave function as a single particle–hole state, yet a superposition of particle-hole states is required to describe the state of the system, the system is characterized by collectivity. We will refer to the electron interaction and the entangled states to which it leads as electron correlation effects. The comparison of the results obtained within the intrachannel and interchannel schemes enables us to distinguish the impact of the configuration interaction, i.e., the electron correlation effects, onto a certain process.

Furthermore, mean-field approaches can be used to simplify the dynamics beyond the intrachannel model. For instance, the Hartree-Fock-Slater (HFS) model Slater (1951) reduces the complexity because the exchange interaction between the electrons is modeled in the local density approximation. This results in an effective one-particle picture where the electron experiences a mean-field potential $V(\mathbf{r})$ created by the other electrons in the shell:

[TABLE]

where $V_{\mathrm{ex}}(\mathbf{r})=-\frac{3}{2}\left[\frac{3}{\pi}\rho(\mathbf{r})\right]^{1/3}$ is the Slater exchange potential. Here $Z$ is the nuclear charge, and $\rho$ is the electron density. For the case of large distances from the origin the Latter tail correction Latter (1955) can be used to obtain the proper asymptotic potential for both occupied and unoccupied orbitals:

[TABLE]

with the effective charge $Z^{\prime}=Z-N_{\mathrm{elec}}$, $N_{\mathrm{elec}}$ being the number of electrons. In this way, we obtain the asymptotic behavior of $Z/r$ for small radii from the ion because the electron can feel the whole nucleus as a Coulomb attractor. Another asymptote, $1/r$, is obtained for large distances, since from afar the $N_{\mathrm{elec}}-1$ positive nuclear charges are screened by the other electrons, such that the electron experiences an effective Coulomb potential equivalent to one unscreened charge. Naturally, this model can only be applied in a meaningful way if electron correlations are expected to be negligible.

III.3 Algorithm for wave-function propagation and for the diagonalization of the Hamiltonian

The Schrödinger equation (5) can be solved numerically by time propagation. Different schemes have been employed for propagating the wave function Feit et al. (1982); Tal-Ezer and Kosloff (1984), but we will focus here on the Lanczos-Arnoldi algorithm Golub and Loan (1996). The well-known fourth-order Runge-Kutta algorithm Hoffman (2001), which is often employed for solving differential equations of first order, is not symmetric in time, but stable if the time step is chosen to be sufficiently small. Since in standard atomic physics problems forward propagation of the Schrödinger equation suffices to obtain the relevant observables and dynamics (e.g. ionization yields, cross sections, electron spectra, etc.) time symmetry is not crucial. If, however, in addition backward propagation is involved in the computations higher accuracy in propagation is needed Goetz et al. (2015). To this end, the Lanczos-Arnoldi algorithm, which is a Krylov subspace method Golub and Loan (1996), can be applied. It is an iterative approximation scheme for eigenvectors and eigenvalues of a matrix, which is, in general, the method of choice for large matrices. It has been used for atomic and molecular physics problems before Kosloff (1994).

III.3.1 Lanczos-Arnoldi algorithm

Suppose that the eigenvalues of a large and sparse Hermitian matrix $A\in\mathbb{C}^{n\times n}$ shall be calculated. Let $\lambda_{1},\ldots,\lambda_{n}$ be the $n$ eigenvalues of $A$, ordered by their magnitude. The Rayleigh quotient

[TABLE]

with $\mathbf{x}\in\mathbb{C}^{n}$, yields the smallest and the largest eigenvalues by the relations $\lambda_{1}=\underset{\mathbf{x}\neq\mathbf{0}}{\max}\,R(\mathbf{x})$ and $\lambda_{n}=\underset{\mathbf{x}\neq\mathbf{0}}{\min}\,R(\mathbf{x})$ Golub and Loan (1996). Let $\mathcal{V}_{k}$ be a subspace of $\mathbb{C}^{n\times n}$ and let $\{{\mathbf{q}}\}_{1}^{k}=\left\{{\mathbf{q}}_{1},{\mathbf{q}}_{2},\ldots,{\mathbf{q}}_{k}\right\}$ be an orthonormal basis of $\mathcal{V}_{k}$. Arranging the vectors ${\mathbf{q}}_{k}$ as columns in a matrix $Q_{k}$ the eigenvalues of $Q_{k}^{T}AQ_{k}$ shall approximate the eigenvalues of $A$. The Lanczos-Arnoldi method generates the vectors ${\mathbf{q}}_{k}$ iteratively, such that the eigenvalues of the matrices $Q_{k}^{T}AQ_{k}=T_{k}\in\mathbb{C}^{k\times k}$, with $k<n$, are, with $k\rightarrow k+1$, progressively better approximations to the eigenvalues of A. Based on the Courant-Fischer theorem Golub and Loan (1996) the algorithm determines increasingly better eigenvalues consider the Rayleigh quotient of the matrix $T_{k}$:

[TABLE]

with $\mathbf{y}\in\mathbb{C}^{n}$ and vectors $\mathbf{u}_{k},\mathbf{v}_{k}\in\mathrm{span}\left\{{\mathbf{q}}_{1},{\mathbf{q}}_{2},\ldots,{\mathbf{q}}_{k}\right\}$. $R(\mathbf{x})$ increases most rapidly in the direction of the gradient

[TABLE]

from which it follows that $\mathbf{\nabla}R(\mathbf{x})\in\mathrm{span}\left\{\mathbf{x},A\mathbf{x}\right\}$. The largest eigenvalue of the next iteration step $M_{k+1}$ will be larger than $M_{k}$, and therefore approach the “real” eigenvalue of the original matrix $A$ if the next vector $\mathbf{q}_{k+1}$ is determined such that $\mathbf{\nabla}R(\mathbf{x})\in\mathrm{span}\left\{{\mathbf{q}}_{1},{\mathbf{q}}_{2},\ldots,{\mathbf{q}}_{k+1}\right\}$. Following the same argument for the eigenvalue of minimal magnitude, if also $\mathbf{\nabla}R(\mathbf{v}_{k})\in\mathrm{span}\left\{{\mathbf{q}}^{1},{\mathbf{q}}^{2},\ldots,{\mathbf{q}}^{k+1}\right\}$ then $m_{k+1}<m_{k}$, because $R(\mathbf{x})$ decreases most rapidly in the direction of the negative gradient $-\mathbf{\nabla}R(\mathbf{x})$. Therefore, both requirements can be satisfied if $\mathbf{q}_{k+1}$ is chosen such that

[TABLE]

and, thereby, successively applying higher powers of the matrix $A$ the Krylov space is built.

In our case we wish to solve the time-dependent Schrödinger equation (5). Formally this equation has the solution

[TABLE]

where

[TABLE]

is the time evolution operator and $\mathcal{T}$ denotes the time-ordering operator

[TABLE]

Approximating the time evolution operator for small time arguments the wave function at the next time step $t+\mathrm{d}t$ has the value

[TABLE]

The Krylov space is built by acting on the starting vector $|\Psi(0)\rangle$ with increasingly higher powers of the Hamiltonian

[TABLE]

where the set of $N$ vectors ${\mathbf{v}}_{n}$, $n=0,1,\ldots N-1$, forms the Krylov basis. $N$ is the dimension of the Krylov space which for numerical reasons should be as small as possible. Now, suppose we had a tridiagonal matrix $T=Q^{T}\hat{H}Q$ with $Q$ being orthogonal, then we could find an orthogonal matrix $U$ ($U^{T}U=\mathbb{1}$) that diagonalizes $T$. Let us call this diagonal matrix $D=U^{T}TU$. Then from $D=U^{T}Q^{T}\hat{H}QU$ it follows that

[TABLE]

and from the properties of diagonal matrices we have

[TABLE]

Orthonormalizing the $N$ vectors $\{{\mathbf{v}}_{n}\}$ we obtain an orthonormal set of vectors $\{\mathbf{q}_{n}\}$. If we arrange them as columns in a matrix $Q$ we know that the matrix $Q^{T}\hat{H}Q=T$ is tridiagonal (QR factorization) Golub and Loan (1996):

[TABLE]

We wish to directly compute the elements of this tridiagonal matrix, $\{\alpha_{k}\}_{k=1}^{N}$, $\{\beta_{k}\}_{k=1}^{N-1}$, in an iterative way. Since $\hat{H}Q=QT$ we find

[TABLE]

while $\beta_{0}{\mathbf{q}}_{0}\equiv 0$. Solving this equation for ${\mathbf{q}}_{k+1}$, if

[TABLE]

then ${\mathbf{q}}_{k+1}={\mathbf{r}}_{k}/\beta_{k}$, where $\beta_{k}=|{\mathbf{r}}_{k}|$. The vectors ${\mathbf{q}}_{k}$ are called Lanczos vectors. In this way the Lanczos-Arnoldi algorithm is obtained, where $\mathbf{q}_{0}$ denotes the starting vector:

[TABLE]

Having calculated the sets of $\alpha$ and $\beta$ coefficients, the tridiagonal matrix $T$ is built up and diagonalized to yield the eigenvalues that form the matrix $D$. Then the matrix product $P\equiv QU$ is calculated. Hence, we can evaluate Eq. (26), the advantage being that we have just to perform matrix-vector multiplications in the small $N$-dimensional Krylov space (typically on the order of $N=10\sim 20$) instead of the whole Hilbert space.

The Lanczos-Arnoldi propagation method yields the same results as the Runge-Kutta algorithm with the additional feature of a significantly higher precision when used in the backward direction Goetz et al. (2015). This algorithm can also be used to calculate approximate eigenstates of the Hamiltonian, which will be exploited in the next section.

III.3.2 Diagonalization of the Hamiltonian: Diabaticity in tunneling ionization

Hamiltonian eigenvalue problems can involve very large matrices. For example, when high energies are involved or strong fields are considered, the necessary number of virtual states can be on the order of $10^{3}$ and the maximum angular momentum can be on the order of $10^{2}$. In combination with many active channels the dimension of the Hilbert space can easily reach $10^{6}$. Typically, these matrices are sparse, but nevertheless an efficient diagonalization algorithm is indispensable. Since the diagonalization of large matrices is one of the most common procedures that are treated numerically it is not surprising that there is software providing a method for the solution of huge eigenvalue problems. The Arnoldi software package ARPACK Lehoucq et al. (1998) employs the iterative Lanczos-Arnoldi algorithm to solve a general eigenvalue problem

[TABLE]

As an application for the diagonalization of the $N$-electron Hamiltonian we will examine the situation when helium is irradiated by an intense electric field $F(t)$ with a small photon energy compared to the electron binding energy. Let the process be characterized by a Keldysh parameter $\gamma\ll 1$, i.e., it is described by tunneling ionization. As mentioned above, for small frequencies the tunneling regime can be viewed as quasi-static, which means, that the system is given time to adjust to the parameters on which it depends (in this case the electric field strength). Such a response is called adiabatic. Therefore, let us analyze the associated adiabatic atomic eigenstates of helium along the lines of Ref. Karamatskou et al. (2013). We would like to compute the eigenstates and the eigenvalues of the Hamiltonian:

[TABLE]

where $E_{n}$ are the eigenenergies corresponding to the eigenvector $|\Psi_{n}\rangle$ of the system. Here, $\hat{H}=\hat{H}_{0}+\hat{H}_{1}$ is the field-free Hamiltonian of Eq. (6), $F(t)$ is the electric field in the $z$ direction and $\hat{z}$ is the position operator. The term $F(t)\,\hat{z}$ describes the dipole interaction between field and electron for a linearly polarized light field in the $z$-direction and thereby represents the light-matter interaction in the length form. The Hamiltonian is equivalent to Eq. (6) as long as no higher multipole orders play a role222Remember, that the vector potential $\mathbf{A}(t)$ and the electric field $\mathbf{F}(t)$ are related: $\mathbf{A}(t)=-\int_{-\infty}^{t}\mathbf{F}(\tau)d\tau$, $\mathbf{F}(t)=-\frac{d}{dt}\mathbf{A}(t).$ Craig and Thirunamachandran (1998). Instead of solving the corresponding time-dependent Schrödinger equation

[TABLE]

the instantaneous, or adiabatic eigenvalue problem can be solved, where for each time $t$ the field $F$ has a certain value. At a given time $t$, the instantaneous eigenstates, which constitute the adiabatic basis, are defined by Eq. (30). To distinguish adiabaticity from the onset of nonadiabatic effects we expand the electronic wave function in the adiabatic basis, $|\Psi(t)\rangle=\sum_{n}\alpha_{n}(t)|\Psi_{n}(t)\rangle$. Upon inserting this expression into Eq. (31) and projecting onto the eigenstate $|\Psi_{m}(t)\rangle$, the equation of motion for the coefficient $\alpha_{m}(t)$ reads

[TABLE]

The off-diagonal matrix elements (also called nonadiabatic coupling terms) $\langle\Psi_{m}(t)|\dot{\Psi}_{n}(t)\rangle$ introduce couplings between different adiabatic eigenstates and lead to nonadiabatic dynamics Zener (1932). If these couplings are very small the so-called adiabatic approximation can be made and Eq. (32) becomes

[TABLE]

This equation shows that in the adiabatic approximation the system remains in a specific adiabatic eigenstate and evolves only with a phase.

Now, it comes handy that for every tunneling state, i.e., every adiabatic state that allows the electron to tunnel through the field-induced barrier, there exists a discrete eigenstate of the instantaneous Hamiltonian (a so-called Siegert state) Bohm et al. (1989); Siegert (1939). Such a state is associated with a complex energy and lies by definition outside of the Hermitian domain of the Hamiltonian because the associated wave function is exponentially divergent for large distances from the atom. One strategy to eliminate the divergent behavior and render the tunneling wave function square integrable is to add a complex absorbing potential (CAP) to the Hamiltonian, which yields complex eigenenergies Riss and Meyer (1993); Santra and Cederbaum (2002). The imaginary part of the Siegert energy $E$ provides the tunneling rate $\Gamma$ of each Siegert state by the relation $\Gamma=-2\ {\mathrm{Im}}(E)$ Moiseyev (1998).

Solving Eq. (30) including a CAP in the Hamiltonian for various field strengths yields the adiabatic eigenstates of the atom and the corresponding (complex) eigenenergies as a function of the field strength. The real part of the Siegert energies are shown in Fig. 4.

Many avoided crossings among the higher adiabatic eigenstates can be observed for field strengths in the range below $0.01$ a.u. ($1$ a.u.$=5.14\times 10^{9}$ V/cm). The ground state energy does not change significantly up to field strengths around $0.02$ a.u. As shown in Ref. Karamatskou et al. (2013), the avoided crossings around $0.02$ a.u. (see inset of Fig. 4) are simply jumped over and the system does not follow the adiabatic ground state (lowest red curve). Rather, for each applied field strength, the electronic state follows the instantaneous eigenstate that has the maximal overlap with the field-free ground state. Therefore, in order to find a more meaningful description of the electronic state a diabatization method is employed whereby the diabatic state $|\Psi_{0}^{(d)}(t)\rangle$ is constructed from the adiabatic basis $\left\{|\Psi_{n}(t)\rangle\right\}$ by requiring a maximal overlap with the field-free state $|\Psi^{(0)}\rangle$. To put it mathematically, for each field strength we seek the adiabatic state $|\Psi_{n}(t)\rangle$, for which $|\langle\Psi_{n}(t)|\Psi^{(0)}\rangle|>|\langle\Psi_{m}(t)|\Psi^{(0)}\rangle|$ for all $m\neq n$. Then, the diabatic state $|\Psi_{0}^{(d)}(t)\rangle$ is set equal to $|\Psi_{n}(t)\rangle$. As can be seen in Fig. 5 c) the diabatic state has a $>90\%$ overlap with the field-free ground state. Figures 5 a) and 5 b) show the real part of the energy and the tunneling rate as a function of the applied electric field. The quadratic behavior of the real part of the energy reflects the Stark effect and allows to read off the polarizability of the system. As expected, the tunneling rate increases considerably for sufficiently high field strengths, which can be understood within the picture visualized in Fig. 1 a): The barrier width decreases because of the stronger field and thereby the tunneling probability is enhanced.

The fact that the system follows one single diabatic state, which essentially maintains the character of the field-free ground state, gives a clearer and more intuitive picture for the explanation of the physics in the tunneling regime which, so far, was studied in the static case.

Furthermore, the diabatic state can be investigated for a dynamic scenario Karamatskou et al. (2013) which involves a pulse with a Gaussian envelope:

[TABLE]

where $F_{0}$ is the peak field strength, $\tau$ is the pulse duration (full width at half maximum, FWHM) and $\omega$ is the field frequency.

We wish to calculate the ionization probability for the diabatic state $|\Psi_{0}^{(d)}\rangle$ when applying this pulse. Moreover, we want to show that this state is sufficient to describe the dynamics in the tunneling regime because the system remains to a large extent in this state. Assuming that we have found the appropriate diabatic basis this particular diabatic state can be described in the diabatic basis by a single coefficient $\alpha_{0}^{(d)}$. Recalling Eq. (33), we can repeat the analysis for the diabatic case, and, neglecting transitions to other diabatic states (i.e., now making a ”diabatic approximation”), we obtain the following equation of motion for the coefficients:

[TABLE]

where $\Gamma_{0}^{(d)}$ is the decay rate of the diabatic state. The population evolution of our distinguished diabatic state, $P_{0}^{(d)}(t)=|\alpha_{0}^{(d)}(t)|^{2}$, during the pulse can be calculated with the equation of motion by taking the time derivative

[TABLE]

where, in the last step, Eq. (35) was inserted. Of course, as we saw above, the ionization rate depends on the external field. This equation can easily be solved analytically by separation of variables:

[TABLE]

with the initial condition $P(t\!=\!-\infty)\!=1$, i.e., all population being in this state long before the pulse. Since Eq. (37) accounts only for the population in our distinguished diabatic state, deviations in the population dynamics can be attributed to nondiabatic behavior, i.e., transitions to other diabatic states. Therefore, we can compare the solution of this equation, in the diabatic approximation giving only the time evolution of our diabatic state, to the solution of the Schrödinger equation, which takes into account all different states and their transitions. The result of this comparison for different photon energies is shown in Fig. 6 for a rather large peak field strength of $0.2$ a.u. Depicted are the populations calculated from Eq. (37) and the Schrödinger equation (31), respectively, [Fig. 6 a)] and the relative difference [Fig. 6 b)] between them after the end of the pulse.

For sufficiently low energies, i.e., deep in the tunneling regime, the total ionization probability is reproduced exactly by considering only the diabatic state (region I). This means that the system remains in this single state and, therefore, it suffices for the description of the dynamics. In the language of the adiabatic representation, the tunneling regime, where the Keldysh parameter $\gamma\ll 1$, is governed by a single diabatic state. For higher photon energies around $1$ eV (region II), the difference increases significantly. In the region where $\gamma\approx 1$ the relative difference of the two methods amounts already to $\approx 10\%$, which is a clear sign of transitions to other diabatic states. This indicates that nondiabatic effects start to become important.

Instead of varying the photon energy for a fixed pulse duration we can also divide the frequency range according to the number of cycles in the pulse. For the highest energies that are presented here multi-cycle pulses are considered [cf. Fig. 7b)]. Going to lower photon energies, around $\approx 0.8$ eV,we reach the few-cycle pulse regime [cf. Fig. 7a)], where nice agreement of the diabatic state ionization and the solution of the Schrödinger equation is obtained. In contrast, tunneling theories such as the ADK theory of tunneling ionization and similar approaches cannot reproduce the correct (diabatic) ionization rate for few-cycle pulses Ammosov et al. (1986).

In order to understand this difficulty we observe that the main difference between few- and multi-cycle fields is the change in field strength between consecutive field oscillations. For a few-cycle pulse the pulse envelope – and with it the ionization rate – changes significantly from one field oscillation to the next. This is depicted in Fig. 7 a). As can be seen in Fig. 7 b), the more cycles there are within the pulse, the less the variation in amplitude between consecutive oscillations. In the framework of ADK theory the ionization rate $\overline{\Gamma}(t)$ is obtained by integrating over one period of the field Bisgaard and Madsen (2004)

[TABLE]

where $\Gamma[F]$ is the instantaneous ionization rate. In the few-cycle limit the ionization rate cannot be simply averaged over one period, because the pulse envelope changes dramatically within one cycle. Of course, for multi-cycle pulses the rate can be averaged over the pulse and is then obtained from Eq. (37) as

[TABLE]

Furthermore, we observe that the ionization probability in region I of Fig. 6 is not constant as a function of photon energy although it is well described by the ionization out of $|\Psi_{0}^{(d)}\rangle$ only. According to our analysis, this pronounced frequency dependence in the few-cycle limit should be attributed to the form of the pulse or, equivalently, to the relation between the number of cycles and the pulse envelope.

In this example we learn by employing the adiabatic representation that the Keldysh parameter is an approximate measure of diabaticity. In the few-cycle tunneling regime a single diabatic state is sufficient to describe the ionization dynamics.

So far, we related the ionization dynamics to the eigenstates of the $N$-electron system. Let us now turn to the wave packet of the outgoing electron and the information it carries about the underlying photoionization process. To this end, a method for the calculation of photoelectron distributions is presented in the following.

IV Calculation of photoelectron distributions within TDCIS

Photoelectron spectroscopy has proven a powerful tool to analyze the processes that happen within complex systems upon irradiation and understand their electronic structure Becker and Shirley (1990). Photoelectron spectra (PES) and photoelectron angular distributions (PAD) contain not only information about the interaction of the electrons with the light field, but also about the electronic correlations in the atomic shell Schwarz (1980). Photoelectron distributions can help to extract information and predict fundamental processes occurring during the interaction with the light pulse, e.g., to steer electrons with light waves Kienberger et al. (2002), which opened the way to attosecond streaking techniques Drescher et al. (2002), to study nonsequential and sequential double ionization in atoms Nikolopoulos et al. (2003), to scrutinize multielectron ionization dynamics Boguslavskiy et al. (2012) and multiphoton excitations of deep shells of atoms Meyer et al. (2010), and to reveal important information about the time-dependence of electron dynamics Krausz and Ivanov (2009); Corkum and Krausz (2007).

IV.0.1 Wave-function splitting method

The calculation of the angular and energy-resolved photoelectron distribution allows for a direct comparison of theoretical predictions with experimental data. In principle, it can be done easily after the pulse is over by projecting the photoelectron wave packet onto the eigenstates of the field-free continuum, which corresponds to a Fourier transform of the wave packet. However, this approach requires large numerical grids in cases where the electrons can travel a long distance while the (possibly long and strong) pulse is still interacting with the system. Therefore, this approach can often not be applied to strong-field problems where the description of the ionized wave packet is particularly challenging due to the nonperturbative interaction between the electrons and the light pulse.

For this reason new schemes, which overcome the obstacle of large grids, were developed for the calculation of photoelectron distributions. One method involves splitting the wave function into an internal and an asymptotic part Tong et al. (2006) where the latter is then analyzed to yield the spectrum. Another method measures the electronic flux through a sphere at a fixed radius which allows to infer the parts of the outgoing wave packet Tao and Scrinzi (2012). Both methods, which have been combined with TDCIS Karamatskou et al. (2014, 2015), yield double-differential photoelectron spectra, i.e., the electron distribution as a function of kinetic energy and ejection angle. In the following the wave-function splitting method within TDCIS is presented briefly Karamatskou et al. (2014).

The wave-function splitting method rests on the assumption that once the ejected electron has traveled far enough from the ion it can be viewed as completely free and merely interacting with the laser pulse that might still be present. This is called the Volkov solution or strong-field approximation Reiss (1980, 2008). Obviously, since in reality the Coulomb potential is a long-range potential the above-mentioned approximation will never be exactly true. Hence, further approximations treating the Coulomb field in higher order, like the Coulomb-Volkov approach were introduced Jain and Tzoar (1978); Duchateau, G. et al. (2000) and studies were performed in order to quantify the influence of the long-range potential on the physical observables Zhang and Nakajima (2007). However, often the potential’s long-range character can be presumed to be negligible (e.g., fast electrons should be influenced only mildly at large distances), which renders this approximation valid.

To calculate the spectral components, the outgoing parts of the wave packet are split off the wave function of the system once they have reached a region far away from the atom. Since during the splitting procedure the laser field is still nonzero these parts must be analyzed using the Volkov Hamiltonian, which describes a free electron in the presence of an electric field Tong et al. (2006).

The theoretical framework of Sec. III is employed to collect for each ionization channel all single excitations from the occupied spin orbital $|\varphi_{i}\rangle$ in one “channel wave function” by taking the sum over all virtual orbitals Rohringer et al. (2006):

[TABLE]

Using these channel wave functions all quantities are calculated in a channel-resolved manner. In this way, effectively, one-particle wave functions for each particular ionization channel $i$ are obtained.

As the method’s name suggests, the channel wave function (40) is split smoothly by applying the radial splitting function Tong et al. (2006)

[TABLE]

The parameter $r_{c}$ determines the center of the splitting function, and $\Delta$ is a “smoothing” parameter controlling the slope of the function. The splitting radius must be chosen sufficiently large, such that the electron is already far enough away to not return to the ion. The splitting function must not alter the ground state, i.e., $\hat{S}|\Phi_{0}\rangle=0$. This is achieved by requiring that $r_{c}/\Delta\gg 1$. During the time propagation of the wave function the splitting is applied at certain splitting times $t_{\mathrm{spl}}$. At the first splitting time $t_{0}$ the channel wave function is split into two parts:

[TABLE]

As visualized in Fig. 8, when the splitting is applied $|\chi_{i,\mathrm{in}}(t_{\mathrm{spl}})\rangle$ is the wave function in the inner region $0<r\lesssim r_{c}$ and $|\chi_{i,\mathrm{out}}(t_{\mathrm{spl}})\rangle$ is the wave function in the outer region $r_{c}\lesssim r\leq r_{\mathrm{max}}$.

We assume that far from the ion the electron experiences only the laser field and not the Coulomb field of the parent ion. Therefore, the outer part of the wave function $|\chi_{i,\mathrm{out}}(t_{\mathrm{spl}})\rangle$ can be propagated analytically to a long time $T$ after the laser pulse is over using the Volkov Hamiltonian:

[TABLE]

which involves the vector potential of the electric field. The corresponding Volkov time propagator has the form

[TABLE]

The inner part of the wave function $|\chi_{i,\mathrm{in}}(t_{\mathrm{spl}})\rangle$ still experiences the Coulomb couplings of the parent ion, and must be propagated using the full CIS Hamiltonian [see Eqs. (11a) and (11b)].

At each splitting time $t_{\mathrm{spl}}$ the outer part of the wave function is extracted. The inner part is propagated with the full Hamiltonian to the next splitting time and subsequently split into inner and outer part, and so forth, until all parts of the ejected electron wave packet have reached the outer region. Each outer part of the wave function is propagated analytically to large times after the pulse.

In order to calculate the spectral components of the emitted electron’s wave packet, $|\chi_{i,\mathrm{out}}(t_{n})\rangle$ is expressed in the CIS basis using new expansion coefficients $\beta_{i}^{a}(t_{n})=\langle\varphi_{a}|\hat{S}|{\chi}_{i}(t_{n})\rangle$ such that wave function in the outer region reads Karamatskou et al. (2014)

[TABLE]

Employing the splitting method the Eqs. (11) have the following form for the inner and outer part of the wave function, respectively Karamatskou et al. (2015):

[TABLE]

The term involving $(2v_{ai^{\prime}ib}-v_{ai^{\prime}bi})$ in Eq. (46a) vanishes for large distances as $1/r$ Rohringer et al. (2006), and, therefore, the Volkov approximation for the outer part of the wave function is justified. The last term in Eq. (46b) couples different channel indices, $i,i^{\prime}$, which introduces channel mixing in the outer part of the wave function. In order to solve this system of differential equations we first observe that the time evolutions of the ionic and the electronic part of the electron wave packet in the outer region are decoupled and can be propagated with two time evolution operators to a large time $T$. Calling the ionic operator $\hat{U}^{\mathrm{ion}}(T,t)$ and the electronic Volkov time propagator $\hat{U}^{\mathrm{elec}}(T,t)=\hat{U}_{V}(T,t)$, see Eq. (44), we can write the time evolution for the coefficients of the outer part of the wave function as Karamatskou et al. (2015)

[TABLE]

Inserting this into Eq. (46b) yields the following equation of motion for the ionic time evolution operator:

[TABLE]

with $U_{ij}^{\mathrm{ion}}(t,t)=\delta_{ij}$. It is solved by numerical time propagation (e.g., using the Lanczos-Arnoldi algorithm) up to large times $T$ and thereby the coefficients $\beta_{i}^{b}(t_{n})$ are obtained.

Now, we need to calculate the spectral components. The Volkov states $|\mathbf{p}\,^{V}\rangle=(2\pi)^{-3/2}e^{i\mathbf{p}\cdot\mathbf{r}}$, which in the velocity form are plane waves, are the eigenstates of the Volkov Hamiltonian and form a basis set in which the channel wave packet for each channel $i$ at final time $T$ can be expanded:

[TABLE]

where a sum of the spectral components over all splitting times $t_{n}$ must be performed. Inserting Eq. (45), the corresponding coefficients are obtained at final time $T$ as

[TABLE]

Then, an incoherent summation over all possible ionization channels must be performed to obtain the double-differential photoelectron distribution as a function of the kinetic energy and the angle with respect to the light polarization axis (because of the linear polarization the distribution is symmetric in the azimuthal angle):

[TABLE]

The extra factor of $p$ results from the conversion from the momentum to the energy differential. The final time $T$ must be chosen sufficiently large, in order to ensure that all parts, in particular the low-energy components, of the outgoing wave packet have reached the outer region and can be analyzed. Further details can be found in Ref. Karamatskou et al. (2014).

V Applications

In the following a few examples of multiphoton ionization of atoms will be presented. In particular we choose the process of above-threshold ionization (ATI) for the investigation in different photon energy regimes. ATI, first observed in 1979 by P. Agostini et al. Agostini et al. (1979), is a highly nonlinear process where an electron absorbs more photons than are necessary for ionization Gontier et al. (1980). As a consequence, a series of peaks can be observed in the photoelectron spectrum, where the separation between two consecutive peaks corresponds to the energy of one photon.

V.1 ATI in the XUV energy range

Let us employ the wave-function splitting method to understand the ATI of argon ($Z=18$, closed shell $1s^{2}2s^{2}2p^{6}3s^{2}3p^{6}$) in the XUV energy range. Numerically, for the splitting method to work properly, three parameters have to adjusted: the splitting radius, $r_{c}$, the smoothness of the splitting function, $\Delta$ [cf. Eq. (41)], and the time instances at which the absorption is applied during the time propagation, $t_{\rm spl}$. In the following we apply a strong Gaussian XUV pulse centered at $\omega=105$ eV with $\Delta\omega=0.7$ eV bandwidth (FWHM), which corresponds to a Fourier-transform limited pulse with $\tau=108$ a.u. ($2.6$ fs) duration (FWHM).333A pulse having the minimum possible pulse duration given an energy bandwidth is called Fourier-transform limited and the relation holds (in atomic units): $\Delta\omega\Delta\tau=2.765$. The peak intensity of the pulse is $1.0\times 10^{15}$ Wcm*-2*. At this photon energy the $3s$ and the $3p$ electrons in argon can be ionized with one photon. If the pulse intensity is high enough they can absorb two photons and undergo ATI.

First we present the double-differential photoelectron distribution, i.e., the ionization yield as a function of energy and angle Karamatskou et al. (2014). In Fig. 9 a) the full angle- and energy-resolved photoelectron distribution of argon after one-photon absorption is shown. The energies of the peaks correspond to the difference between the photon energy and the binding energy of the corresponding orbital ($3s$ and $3p$, respectively). The energy width of the peaks corresponds to the Fourier-transform limited energy width. Fig. 9 b) shows the corresponding ATI peaks, which are separated in energy from the one-photon peaks exactly by the energy of one photon.

The second axis represents the angle with respect to the polarization axis. The angular distributions, which are different for the different peaks, feature the corresponding contributions from the different channels. They are visualized in Fig. 9 c) as cuts along the maxima of the corresponding peaks. In one-photon ionization the $3s$ electron occupies a $p$ state ($l=1$) which manifests itself in the $Y_{10}\propto\cos^{2}\theta$ distribution. The $3p$ one-photon peak has both an $s$- and a $d$-wave contribution. Analogously, through the absorption of two photons the ATI peak of the $3s$ electron exhibits an $s$- and $d$-wave character and the $3p$ peak a $p$- and $f$-wave character.

When integrated over the solid angle the photoelectron distribution is just a function of the photoelectron’s kinetic energy:

[TABLE]

In that case it is commonly called photoelectron spectrum (PES). In Fig. 10 the photoelectron spectrum corresponding to the photoelectron distribution of Fig. 9 is shown, i.e., integrated over the solid angle. The peaks reflect the binding energy position and allow for the characterization of the corresponding shell that was ionized. Such energy-dependent PES are recorded in experiment when all electrons irrespective of their emission direction are measured, e.g., by a magnetic-bottle spectrometer. In the next subsection, this type of spectra measured for atomic xenon will be of central importance.

V.2 Xenon ATI involving the giant dipole resonance

Now we will investigate a system exhibiting electron correlations in the atomic shell by employing multiphoton ionization as a tool.

In 1964 a strong response of xenon ($Z=54$) to XUV radiation was discovered in the one-photon absorption spectrum Ederer (1964), the so-called giant dipole resonance (GDR). The three outermost shells of xenon are the $5p^{6}$, $5s^{2}$, and $4d^{10}$ shells. The GDR was soon attributed to the interplay of two effects concerning the $4d$ shell. Qualitatively the phenomenon can be explained as originating from a shape resonance effect: After absorbing one photon a $4d$ electron is promoted to the continuum. Due to the centrifugal barrier (predominantly angular momentum $l=3$), it is trapped temporarily in a resonance state near the ionic core Cooper (1964); Starace (1982) until it tunnels out and leaves the ion. The XUV resonance is interpreted as the collective response of all ten $4d$ electrons Fetter and Walecka (1980). Only when electron correlation effects within the $4d$ shell are included is quantitative agreement with experimental data achieved. This means that the resonant excitation cannot be explained as a purely independent-particle effect. Until recently it was assumed that the GDR consists of one single, broad resonance which is associated with one quantum state. In this section it will be shown that using XUV two-photon spectroscopy, there is substantial sensitivity to substructure of this resonance. In fact, the GDR accommodates two resonances.

The GDR has been studied previously within TDCIS Pabst and Santra (2013); Krebs et al. (2014). In Fig. 11 the one-photon cross section of xenon is shown for the intrachannel and interchannel models. The interchannel model captures the main features of the many-body effect that renders the one-photon cross section curve broader and shifts it to higher energies. In fact, the interchannel curve reproduces reasonably well both the position and the width of the experimental cross section Krebs et al. (2014).

V.2.1 Nonlinear response regime

Let us now review the analysis of the GDR in the two-photon regime and the comparison of the theoretical results with experimental data obtained at the free-electron laser FLASH, which was originally published in Ref. Mazza et al. (2015). In the experiment xenon gas was irradiated by an FEL beam and the resulting photoelectrons were recorded with a magnetic-bottle electron spectrometer. As we saw above, photoelectron spectra allow to disentangle one-photon ionization from ATI. One of the spectra that were recorded is shown in Fig. 12.

The processes of interest are the one- and two-photon ionization of the $4d$ shell. Since the photon-energies selected in the experiment lie exactly in the range of the GDR, $105$ eV and $140$ eV, the two-photon process occurs through the giant dipole resonance as an intermediate step (Fig. 13). Subsequently the inner-shell vacancy decays via Auger decay Werme et al. (1972). In this case the two-photon process is ATI because the photon energies exceed the binding energy of the $4d$ orbital. In Fig. 13 the levels of xenon and the ionization processes are shown schematically. All of these processes can be identified in the experimental spectra, cf. Fig. 12. At kinetic energies around $140$ eV, the two-photon ionization from the $4d$ shell is observed in a spectral feature which resembles the shape of the $4d$ (one-photon) emission lines around $35$ eV and is separated from them by exactly the energy of one photon.

In order to quantify the ionization processes the two-photon cross section is calculated within TDCIS by calculating the depopulations in the $4d$ shell due to one- and two-photon absorption Greenman et al. (2010). As long as perturbation theory is valid and higher order processes (e.g., three-photon processes) are negligible, the depopulations are distinguishable due to the different angular momenta of the ejected electron according to the dipole selection rules. All contributions from the $4d_{m}$ subshells, with $m$ being the magnetic quantum number, i.e., $4d_{0}$, $4d_{\pm 1}$, and $4d_{\pm 2}$, must be added. Equivalently, the corresponding peak of the photoelectron distribution can be integrated over the angle and the energy in order to obtain the ionized population.

In the following, the role of collectivity shall be investigated in the nonlinear regime. As described in Sec. III a comparison between the interchannel and intrachannel models can elucidate the impact of collectivity. For the description of a collective response the system cannot be written as a single particle-hole state, but instead a superposition of particle-hole states is needed. The interchannel model includes the coupling among the holes in the $4d$, $5s$ and $5p$ orbitals and the electron, see Fig. 13 d), cf. also Fig. 3 b). The corresponding Coulomb matrix elements $\langle\Phi_{i}^{a}|\hat{V}_{\mathrm{e-e}}|\Phi_{j}^{b}\rangle$ between the particle–hole excitations $|\Phi_{i}^{a}\rangle,\,|\Phi_{j}^{b}\rangle$, are included for all different hole pairs $(i,j)$ within the space of active orbitals ($4d$, $5s$, and $5p$), which means that $\langle\Phi_{i}^{a}|\hat{V}_{\mathrm{e-e}}|\Phi_{j}^{b}\rangle\neq 0$, for all $i=j$ and $i\neq j$; $a,b$ are taking values for all virtual orbitals. In this way, superpositions of particle-hole states, i.e. collective states, may be described. In contrast, in the case of the intrachannel model the elements with $i\neq j$ are set to zero. Keeping only the elements $\langle\Phi_{4d_{m}}^{a}|\hat{V}_{\mathrm{e-e}}|\Phi_{4d_{m}}^{b}\rangle\neq 0$ results, therefore, in the description of coupling only with the $4d_{m}$ orbital from which the electron was ionized, see Figs. 13 c), cf. also 3 a). The relevant cross sections are calculated for both schemes and in the following the results will be compared to the experimental data.

As introduced in Sec. III rate equations can describe the ionization probability in the perturbative regime. The corresponding differential equations describing the ionization out of the $4d$ shell by a light pulse have the form:

[TABLE]

where $P_{0}$ denotes the $4d$ ground state population. The first equation describes the change in population of the ground state due to one- and two-photon ionization, denoted by $p_{1ph}$ and $p_{2ph}$, respectively. $\sigma^{(1)}$ and $\sigma^{(2)}$ are the one-photon and ATI cross sections, respectively. $j(t)=I(t)/\omega$ is the time-dependent flux, where $I(t)$ is the pulse intensity and $\omega$ is the photon energy. In other words, the flux is a measure of the number of photons that are available for ionization processes per unit time and area.

The rate equation system (53) is solved for a variety of pulse intensities spanning many orders of magnitude (e.g., using the Runge-Kutta algorithm). One- and two-photon ionization yields are obtained for both the intrachannel and interchannel model by inserting the corresponding cross sections.

In our considerations the ground state population can be set constant, i.e., $P_{0}(t)=P_{0}=1$ because a crucial criterion for the perturbative regime is that the ground state is not depleted. Thereby the equations are simplified and can be integrated to yield the depopulations due to $N$-photon absorption

[TABLE]

Here, $F^{(N)}=\int j^{N}(t)dt$ denotes the fluence that is available for the ionization by $N$ photons. For all our purposes we will use Gaussian pulses with the intensity envelope

[TABLE]

where $c$ is the speed of light, $E_{0}$ is the peak electric field, and $\tau$ is the pulse duration (FWHM). We obtain the fluences for one-photon and two-photon absorption:

[TABLE]

From Eq. (54) the expressions for the cross sections follow: $\sigma^{(N)}=P_{N}/F^{(N)}$. The units of (generalized) cross sections are cm2NsN-1 Saenz and Lambropoulos (1999).

V.2.2 Comparison between theory and experiment

In order to meaningfully compare the theoretical calculations with experimental data it is necessary to perform a volume integration matching the experimental setting. Because of the experimental spatial pulse profile the light intensity varies over the interaction region of the light with the gas. Therefore, also the ionization probability depends on the position of an atom within the electric field distribution. To find the correct spatial dependence the following parameters are important: the duration and the statistics of the light pulses, the beam focus, the Rayleigh length, the volume of acceptance along the propagation direction and perpendicular to it, described by the coordinate $z$ and the radial coordinate $\rho$, respectively.

To this end, we use our signal, $S(F)$, which is given as a function of intensity or, equivalently, of fluence. A Gaussian beam profile is assumed, i.e., the fluence distribution is given by $F(\rho,z)=F_{0}(z)\exp[-\rho^{2}/w(z)^{2}]$, with $F_{0}(z)=4n_{phot}\ln 2/[\pi w^{2}(z)]$. $n_{phot}$ is the number of photons in the pulse. The function $w(z)$ describes the divergence of the beam, i.e. it determines the spot size at a distance $z$ from the beam waist or focus $w_{0}$:

[TABLE]

Denoting the integration volume element by $\mathrm{d}V$ we obtain the expression

[TABLE]

where a coordinate transformation $(\rho,z)\rightarrow(F,z)$ is performed with the determinant of the Jacobian transformation matrix

[TABLE]

because an integration over the fluence is more favorable than over the radial coordinate. The fluence distribution can be inverted bijectively and yields $\rho(F,z)$ when we take the positive (physical) solution. We finally obtain

[TABLE]

which can be calculated numerically employing the experimental parameters.

In Fig. 14 the comparison between the two theoretical models, intrachannel and interchannel, and the experimental data is drawn for the two photon energies of $105$ and $140$ eV Mazza et al. (2015). It clearly shows that the interchannel model (solid lines) reproduces the intensity dependence of the experimental yields, whereas the intrachannel model (dashed lines) fails to do so. Note that the slopes of the dash-dotted lines represent the intensity dependence of a linear and a quadratic process according to Eq. (4) and match the experimental points of the corresponding ionization order. The excellent agreement in the ratio between the two orders of ionization and the onset of saturation confirms that the interchannel model captures the relevant physics over the whole intensity range. This leads to the conclusion that correlations between all possible electron-hole states play a major role in the two-photon ionization process. At $140$ eV the experimental results are described particularly well by the interchannel model which gives a much larger cross section than the intrachannel model.

V.2.3 Theoretical analysis of the ATI cross section

As a next step, the influence of collectivity on the cross sections is investigated over the whole photon-energy range of the GDR. In Fig. 15 the one-photon and generalized two-photon cross sections are plotted together for the intrachannel (dashed lines) and the interchannel model (solid lines). The two-photon cross section for the intrachannel model is red-shifted and narrower than the one-photon curve which is indicated by the red arrow spanning the FWHM. This would also be a naive guess for the two-photon cross section, if we consider that the second photon induces a continuum-continuum transition. To understand this, let us assume that a single intermediate state, namely the GDR, is populated resonantly in the two-step ionization process, see Fig. 13 b). In perturbation theory, where the interaction of the atom with the light field, $\hat{H}_{\mathrm{int}}$, is treated as a perturbation, $\sigma^{(2)}$ can be factorized into two one-photon cross sections as long as the photon energy lies is the vicinity of a single, isolated one-photon resonance. The first photon excites the GDR, the second photon ionizes the electron. Recalling Eq. (1), we obtain the two-photon cross section with the GDR as one single intermediate state

[TABLE]

Now, the second photon initiates a continuum-continuum transition. The cross section of this transition follows a simple energy dependence of $E^{-l-7/2}$, cf. Sec. $70$ of Ref. Bethe and Salpeter (1977), with $l$ being the angular momentum of the initial state. The GDR exhibits mostly $f$-character, such that the exponent becomes $-13/2$. Therefore, for one single intermediate state the two-photon cross section is obtained by multiplying the one-photon ionization cross section with this energy-dependent factor. According to the two-step picture one expects a narrower two-photon peak that is shifted to lower energy, because the one-photon cross section for exciting an electron from the intermediate state into the continuum decreases monotonically with increasing energy. This is visualized in Fig. 16 a). The intrachannel one-photon ionization cross section (black curve) multiplied by $E^{-13/2}$ is indeed shifted to smaller energies and exhibits also a smaller width, just like the two-photon cross section peak (red curve). Therefore, the behavior of the two-photon cross section in the intrachannel case can be qualitatively understood in terms of a sequential process involving the GDR as a single intermediate state (dashed blue curve).

Coming back to Fig. 15 we observe a fundamentally different behavior for the interchannel case: Unexpectedly, it predicts a significantly broader two-photon cross section curve (solid red curve) than for the one-photon case, demonstrated by the red arrow which spans only half the FWHM. Moreover, the shape of the curve is peculiar, it does not resemble a simple Lorentzian profile, but exhibits a knee as can be seen more clearly in the inset of the figure. What is the origin of the broadening and the shape?

To answer this question let us apply the two-step model to the interchannel scenario. As shown in Fig. 16 b) the resulting curve for the generalized two-photon absorption cross section within the simplified model (blue curve) is indeed shifted to a smaller energy. But the width is also decreased in strong contrast to the calculated two-photon cross section curve in the full model, which leads to the broadened red curve. Furthermore, the dashed blue curve underestimates the experimental cross section especially at $140$ eV by a considerable factor. This demonstrates that the simple model does not capture the physics of the full model if only a single resonance is taken into account as the intermediate state.

Besides its broad profile, the shape of the ATI curve strongly suggests that there are two resonances underlying the broad GDR, one being located at the peak of the curve and the other one near the knee-type structure. In the case of two resonances there must be a sum over the corresponding resonance states in the cross section expression (61)

[TABLE]

From this expression it becomes clear that once there are more than one intermediate state, interference terms between overlapping resonances arise, whose relative phase can result in a broadening and can change the shape of the cross section curve. Clearly, due to the cross terms in the ATI cross section, the nonlinear process represents a more sensitive observable for testing the hypothesis of two resonances Mazza et al. (2015).

The ATI cross section provides the first hint to the fact that there exist two underlying resonances. In order to investigate the resonance states further, the same strategy as in Sec. III.3.2 is employed: In Ref. Chen et al. (2015) the resonance states, which are similar to the tunneling states we encountered previously, are investigated in detail by diagonalizing the $N$-electron Hamiltonian. Similarly to our considerations above, the resonance states are exponentially divergent in the asymptotic region Siegert (1939). Hence, it is difficult to access them by a Hermitian Hamiltonian. Again, in order to overcome this obstacle the Hamiltonian is rendered non-Hermitian; this time by the exterior complex scaling method (ECS) Moiseyev (1998) (instead of a CAP as done above). Thereby the associated resonance wave functions are transformed into square-integrable functions.

The states of interest are accessible from the ground state by the absorption of one photon. Therefore, the relevant one-photon resonances are found by imposing the requirement that they possess a large overlap with the ground state coupled by one dipole step Chen et al. (2015).

When doing so, three resonances are found in the intrachannel case, which form a group, cf. Tab. 1. They correspond to the three $4d_{\pm m}$, $|m|=0,1,2$, channels. The real part of the resonance group, $\approx 77$ eV, is consistent with the peak position and enhanced magnitude in the one-photon cross section calculated within the intrachannel model, cf. the black-dashed curve in Fig. 15. Activating the interchannel coupling, however, reveals two distinct underlying resonance states as shown in Tab. 1 Chen et al. (2015). This means, that many-body correlations lead to resonances that cannot be attributed to single ionization pathways. The positions of the resonances, marked by arrows in the inset of Fig. 16 at $\approx 74$ and $108$ eV, are also consistent with our findings concerning the ATI cross section. The wave function of the excited electron shows a prominent $f$-wave character for both resonances such that the xenon GDR is dominated by $4d\rightarrow\epsilon f$ transitions. This legitimates our previous assumptions for the energy dependence in the two-step ionization model for the two-photon cross sections. More details can be found in Ref. Chen et al. (2015).

Summarizing, this example shows that the nonlinear response of a many-electron system to intense XUV radiation can be used to unveil information about its collective behavior. The theoretical xenon ATI cross section exhibits a knee-type structure which is not visible in the one-photon cross section. The agreement of the TDCIS model with experimental results in the two-photon regime legitimates the prediction of two resonance states underlying the GDR.

V.3 ATI in the x-ray regime

We continue our investigation towards smaller wavelengths in the x-ray regime. Following the objective of Ref. Tilley et al. (2015) we ask how important ATI is at hard x-ray photon energies. As described above, PES are most adequate observables for quantifying the impact of ATI.

X-ray free-electron lasers (XFELs) provide ultrashort (hard) x-ray pulses at very high intensities. Due to their small wavelength x-rays provide high spatial resolution down to a few Ångströms. Such pulses are of particular interest for the purposes of single-molecule imaging via coherent x-ray scattering at atomic resolution Chapman (2009); Neutze et al. (2000) but also for the investigation of electronic dynamics in atoms and (bio-) molecules on a time scale between attoseconds and tens of femtoseconds Giles (2002); Krausz and Ivanov (2009).

Even though the interaction probability of x-rays with matter is low Thompson and Vaughan (2001), in the high-intensity regime it might be necessary to consider the importance of nonlinear processes affecting electronic dynamics of atomic, molecular or solid-state target systems. With higher photon energy the probability for photoabsorption by electrons of the valence shell decreases significantly. In the x-ray regime the absorption probability for electrons in the valence shells at low intensities is negligible compared to visible or XUV light. For instance, already at $1$ keV a $2p$ electron in a carbon atom absorbs with a cross section of only $10^{-4}$ Mb. This cross section is $3$ orders of magnitude smaller than for the core electrons Yeh and Lindau (1985). The core electrons are more likely to absorb because they have a larger binding energy. When it comes to very high intensities inner-shell electrons might absorb even more than one photon despite of the high photon energy. Especially for imaging experiments at XFELs this circumstance could pose a problem because after diffraction from a sample the real space image needs to be retrieved from the image in the momentum-transfer space. For brighter illumination and a higher signal intensity higher pulse intensities are needed, e.g., in order to image the interior of a virus. One strategy to achieve significantly higher intensities at XFELs is to focus the pulses down to a few nanometers. ATI at such high intensities might be a process that can produce considerable signal in the photon energy regime used in imaging experiments.

Therefore, let us examine the role and the magnitude of nonlinear effects in the x-ray regime under high-intensity conditions that might become available soon at XFELs. The atomic potential is treated on the HFS level, which, as discussed in Sec. III, significantly reduces the computational cost as it spares the calculation of the exact Coulomb interaction between the electrons. The comparison with full TDCIS calculations suggests that, indeed, in this photon-energy regime electron correlations are of minor importance Tilley et al. (2015). Nevertheless, the investigation of the interaction of x-rays with atoms is computationally challenging. The grid size, the parameters for the calculation of the PES using the wave-function splitting method and the propagation time step have to be adjusted in order to numerically capture the high-energy components, i.e., the fast oscillations, of the wave packet in time and space Tilley et al. (2015).

The applicability of the PES calculations within TDCIS in the x-ray regime is demonstrated for hydrogen by comparison with previous work Florescu et al. (2011), which shows nice agreement. Although under the conditions of short-wavelength x-rays our assumption of the dipole approximation (cf. Sec. III) may no longer be valid, it is found that the results underestimate the nondipole results merely by a factor of $2-3$ Tilley et al. (2015). Nevertheless, going beyond the dipole approximation in the light-matter interaction would be highly desirable and remains an interesting topic of future investigations. Furthermore, it is found that it is appropriate to neglect the relativistic correction for the ionized electron.

Of course, as mentioned several times in this tutorial, in order to calculate a meaningful cross section we must make sure to be in the perturbative limit, i.e., the ionized populations remain very small. A Gaussian pulse of $10^{20}$ Wcm*-2* intensity and $0.12$ fs pulse duration is assumed at three typical hard x-ray photon energies of $8$, $10$, and $12$ keV. The cross sections are calculated by integrating the corresponding photoelectron peaks and dividing by the fluence, according to Eqs. (56). The two-photon $K$-shell ATI cross sections for hydrogen and for the chemical elements that are commonly found in organic molecules, namely, carbon, nitrogen, and oxygen are presented in Table 2 Tilley et al. (2015). We expect an increase in the ATI cross section for higher $Z$ values because the electron binding energies are larger. This is indeed the case. Also, the cross section should decrease with larger photon energy because the ratio between the electron’s binding energy and the photon energy decreases, and the electron appears almost to be free.444Recall that it is impossible for a free electron to absorb a photon, because energy and momentum cannot be conserved simultaneously.

Fig. 17 shows the PES of carbon at a photon energy of $10$ keV and a rather high intensity of $10^{24}$ Wcm*-2* for the angles $\theta=0$, $\pi/6$, $\pi/3$, and $\pi/2$ Tilley et al. (2015). The one-photon ionization peak and the first two ATI peaks are shown. Each peak consists of two subpeaks, the one at lower energy being associated with ionization from the $K$ shell, the other one with $L$-shell ionization. At this intensity the depopulation due to $K$-shell ATI in the direction $\theta=0$ becomes comparable and even higher than valence one-photon ionization.

Next, the ionization probability of carbon for the cases of one-photon ionization out of the $L$ shell and of two-photon ionization out of the $K$ shell are presented for different pulse intensities. They are obtained by integrating the corresponding peaks over the energy and the angle as described in Ref. Tilley et al. (2015). In Fig. 18 the ionization probabilities are shown as a function of intensity on a double-logarithmic scale. We recognize the characteristic quadratic behavior of the two-photon ATI peak and the linear behavior in the one-photon valence ionization probability from the slopes of the lines, remember that $\ln P_{N}\propto N\ln I$ from Eq. (4). Because of the high photon energy, deviations from the perturbative ionization do not play a role for the intensity regime below $10^{24}$ Wcm*-2*. As mentioned above, between $10^{23}-10^{24}$ Wcm*-2* the ionization probability due to $K$-shell ATI peak is comparable to the probability to ionize carbon with one photon out of the valence shells.

Re-examining Fig. 17, a small peak around zero kinetic energy is observed, e.g., in the direction $\theta=0$. What is the origin of these slow electrons? As mentioned above, for the splitting method the final propagation time $T$ must be large enough to detect all slow electrons in the splitting region. Since the focus of the study lay on the fast electrons the propagation time was chosen shorter than necessary to collect all slow electrons. Therefore, the height of the peak shown in Fig. 17 probably underestimates the electron yield. A new calculation for a pulse at $10$ keV photon energy, $12$ as pulse duration, and $3.5\cdot 10^{22}$ Wcm*-2* intensity allows to observe the slow-electron peak in its full height. The resulting PES are shown in Fig. 19 for different angles. The height of the slow-electron peak is comparable to the two-photon ATI peak in the direction of $\pi/2$. This suggests that the two processes are of the same, namely, of second order. We can understand this by realizing that the pulse, being short in the time domain, has a large bandwidth in the energy domain which leads to the energy width of the photoelectron peaks. For a pulse duration of $12$ as the energy bandwidth amounts to $\approx 210$ eV for the two-photon peak. When a valence electron absorbs one photon at a certain energy and re-emits a photon of slightly lower energy within the pulse bandwidth, slow electrons carrying away the excess energy are produced. Of course, the absorption and subsequent emission of one photon is also a second order process.

We conclude from this example that ATI remains negligible for intensities at the most recent XFEL experiments with hard x-rays. However, with photon energies at around $10$ keV or below, when entering the regime around $10^{23}$ Wcm*-2* the two-photon ionization probability of the core electrons reaches the same order of magnitude as one-photon valence ionization.

VI Conclusion

We studied the ionization of many-electron systems in various photon-energy regimes, spanning the infrared up to the hard x-ray regime within the TDCIS scheme, focussing on the nonlinear response regime. Light-atom interactions at intensities high enough to induce multiphoton absorption processes were investigated and the information that photoelectron distributions provide was examined. In particular, xenon was studied as a prime example of an atomic system exhibiting strong correlation effects in the atomic shell. Using xenon as a model system it was shown that nonlinear spectroscopy can be employed as a tool to broaden our knowledge about the collective resonance behavior of atomic systems. With the help of TDCIS, which captures electron correlation effects in the atomic shell, the substructure in the GDR of xenon arising from collective effects in the $4d$ shell was uncovered. The underlying resonance states were characterized in detail by analyzing the eigenstates of the Hamiltonian. Furthermore, the diagonalization of the strong-field Hamiltonian in the tunneling regime revealed, that, strictly speaking, tunneling ionization is a nonadiabatic process and that in the case of few-cycle pulses the dynamics can be fully described by one single diabatic state. It was shown that ATI in the x-ray regime does not play a crucial role at current FEL intensities. However, when increasing the pulse fluence further nonlinear ionization might start to play a role.

Summarizing, the TDCIS scheme has found many applications for the study of multiphoton ionization in the tunneling and perturbative multiphoton regimes employing also the calculation of photoelectron distributions. Particularly for the analysis of electron correlation effects it has proven beneficial and the applicability of the methods presented in this tutorial is by far not exhausted.

VII Acknowledgment

This work was supported by the Hamburg Centre of Ultrafast Imaging through the Louise-Johnson Fellowship.