QPROP with faster calculation of photoelectron spectra

TL;DR

This paper introduces QPROP 3.0, an improved computational tool for calculating photoelectron spectra more efficiently using the i-SURFV method, reducing post-processing time in strong-field laser-atom simulations.

Contribution

QPROP 3.0 incorporates the novel i-SURFV method for single-step post-propagation, enhancing the efficiency of photoelectron spectra calculations in strong-field physics.

Findings

i-SURFV reduces post-propagation computational time.

Enhanced parameter handling improves user experience.

Examples demonstrate the efficiency of the new method.

Abstract

The calculation of accurate photoelectron spectra (PES) for strong-field laser-atom experiments is a demanding computational task, even in single-active-electron approximation. The QPROP code, published in 2006, has been extended in 2016 in order to provide the possibility to calculate PES using the so-called t-SURFF approach [L. Tao, A. Scrinzi, New J. Phys. 14, 013021 (2012)]. In t-SURFF, the flux through a surface while the laser is on is monitored. Calculating PES from this flux through a surface enclosing a relatively small computational grid is much more efficient than calculating it from the widely spread wavefunction at the end of the laser pulse on a much larger grid. However, the smaller the minimum photoelectron energy of interest is, the more post-propagation after the actual laser pulse is necessary. This drawback of t-SURFF has been overcome by Morales et al. [F. Morales, T. Bredtmann, S. Patchkovskii, J. Phys. B: At. Mol. Opt. Phys. 49, 245001 (2016)] by noticing that the propagation of the wavefunction from the end of the laser pulse to infinity can be performed very efficiently in a single step. In this work, we introduce QPROP 3.0, in which this single-step post-propagation (dubbed i-SURFV) is added. Examples, illustrating the new feature, are discussed. A few other improvements, concerning mainly the parameter files, are also explained.

Figures

Tables

Table 1. Table 1: Programs names in Qprop  2.0 and 3.0.

Qprop  2.0	Qprop  3.0	Task
hydrogen_im	imag_prop	Imaginary-time propagation
hydrogen_re	real_prop	Real-time propagation
	isurfv	${\left({\hat{H}}_{\ell }-E\right)}^{-1}$ to ${\varphi }_{\ell m}$, ${\partial }_r{\varphi }_{\ell m}$
eval-tsurff(-mpi)	tsurff(_mpi)	Spectrum with t-SURFF

Table 2. Table 2: Matching of the parameter names in the propagate.param file in the vortex and lissajous examples to the variables in ( 4.2 ).

Parameter	Eq. (4.2)	Domain	Default in
			vortex	lissajous
omega-1	${\omega }_1$	${ℝ}_{+}$	$0.114$	$0.114$
omega-2	${\omega }_2$	${ℝ}_{+}$	$0.190$	$0.228$
max-electric-field-1-x	$E_{1x}$	${ℝ}_{0+}$	$0.0267$	$0.0534$
max-electric-field-2-x	$E_{2x}$	${ℝ}_{0+}$	$0.0267$	$0.0$
max-electric-field-1-y	$E_{1y}$	${ℝ}_{0+}$	$0.0267$	$0.0$
max-electric-field-2-y	$E_{2y}$	${ℝ}_{0+}$	$0.0267$	$0.0534$
num-cycles-1	$n_{c1}$	${ℝ}_{0+}$	$6.0$	$8.0$
num-cycles-2	$n_{c2}$	${ℝ}_{0+}$	$10.0$	$16.0$
phase-1-x	${\varphi }_{1x}$	$ℝ$	$0.0$	$0.0$
phase-2-x	${\varphi }_{2x}$	$ℝ$	$0.0$	$0.0$
phase-1-y	${\varphi }_{1y}$	$ℝ$	$0.5\pi$	$0.0$
phase-2-y	${\varphi }_{2y}$	$ℝ$	$1.5\pi$	$0.0$
tau-delay	$\tau$	${ℝ}_{0+}$	$331.0$	$0.0$

Equations

$$\mathrm{i}{\partial }_t|\Psi (t)⟩=\hat{H}(t)|\Psi (t)⟩$$

$$\hat{H}(t)=-\frac{\Delta }{2}-\mathrm{i}\mathbf{A}(t)\cdot \nabla +U(r)-\mathrm{i}{V}_{\text{ im}}(r).$$

$$⟨\mathbf{r}|\Psi (t)⟩=\frac{1}{r}\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{\varphi }_{\ell m}(r,t)Y_{\ell m}(\Omega ).$$

$$⟨\mathbf{r}|\Psi (t)⟩=\frac{1}{r}\sum _{\ell =0}^{N_{\ell }-1}\sum _{m=-\ell }^{\ell }{\varphi }_{\ell m}(r,t)Y_{\ell m}(\Omega ).$$

$$⟨\mathbf{r}|{\Psi }_{m_0}(t)⟩=\frac{1}{r}\sum _{\ell =|m_0|}^{N_{\ell }-1}{\varphi }_{\ell m_0}(r,t)Y_{\ell m_0}(\Omega ).$$

$$V_{\mathrm{im}}=V_{\mathrm{im},\max }{\left(\frac{r-R_{\mathrm{im}}}{W_{\mathrm{im}}}\right)}^{16}\Theta (r-R_{\mathrm{im}}),$$

$$U(r)=\{\begin{array}{ll}-1/r& \text{ if }r<R_{\mathrm{co}}\\ -(2R_{\mathrm{co}}-r)/R_{\mathrm{co}}^2& \text{ if }{R}_{\mathrm{co}}<r<2R_{\mathrm{co}}\\ 0& \text{ if }r>2R_{\mathrm{co}}\end{array}.$$

$$a_{\text{I}}(\mathbf{k},T)=⟨\mathbf{k}(T)|\Theta (r-R_{\mathrm{I}})|\Psi (T)⟩.$$

$$a_{\text{I}}(\mathbf{k},T)={\int }_0^{T}dt{\partial }_t⟨\mathbf{k}(t)|\Theta (r-R_{\mathrm{I}})|\Psi (t)⟩+⟨\mathbf{k}(0)|\Theta (r-R_{\mathrm{I}})|\Psi (0)⟩.$$

$$a_{\text{I}}(\mathbf{k},T)=\mathrm{i}{\int }_0^{T}dt⟨\mathbf{k}(t)|\left[\hat{H},\Theta (r-R_{\mathrm{I}})\right]|\Psi (t)⟩.$$

$$a_{\text{I}}(\mathbf{k},T)$$

$$\delta a_{\text{I}}(\mathbf{k},{\tau }_p,T)$$

$$\displaystyle\delta a_{\textmd{\tiny{I}}}({\bf k},\tau_{p},T)=-\mathrm{i}\int_{\tau_{p}}^{T}dt~{}\int d\Omega~{}\int drr^{2}~{}\Big{(}\delta(r-R_{\mathrm{I}})\psi^{*}_{\bf k}({\bf r},t)\partial_{r}\Psi({\bf r},t)$$

$$\displaystyle+\psi^{*}_{\bf k}({\bf r},t)\Psi({\bf r},t)\frac{\Delta}{2}\Theta(r-R_{\mathrm{I}})\Big{)}$$

$$\delta a_{\textmd{\tiny{I}}}({\bf k},\tau_{p},T)=-\frac{\mathrm{i}R_{\mathrm{I}}^{2}}{2}\int_{\tau_{p}}^{T}dt~{}\int d\Omega~{}\Big{(}\partial_{r}\psi^{*}_{\bf k}({\bf r},t)\Psi({\bf r},t)-\psi^{*}_{\bf k}({\bf r},t)\partial_{r}\Psi({\bf r},t)\Big{)}\Big{|}_{r=R_{\mathrm{I}}}.$$

$$⟨\mathbf{k}(t)|\mathbf{r}⟩={\psi }_{\mathbf{k}}^{\ast }(\mathbf{r},t)=\sqrt{\frac{2}{\pi }}{e}^{\mathrm{i}{k}^{2}t/2}\sum _{\ell m}(-\mathrm{i}{)}^{\ell }{j}_{\ell }(kr)Y_{\ell m}^{\ast }(\Omega )Y_{\ell m}({\Omega }_k)$$

$$\delta a_{\text{I}}(\mathbf{k},{\tau }_p,T)=\sum _{\ell m}\delta a_{\text{I},\ell m}(k,{\tau }_p,T)Y_{\ell m}({\Omega }_k),$$

$$\delta a_{\text{I},\ell m}(k,{\tau }_p,T)$$

$$\int_{\tau_{p}}^{\infty}dt~{}e^{\mathrm{i}Et}\phi_{\ell m}(r,t)=\int_{\tau_{p}}^{\infty}dt~{}e^{\mathrm{i}Et-\mathrm{i}\hat{H}_{\ell}(t-\tau_{p})}\phi_{\ell m}(r,\tau_{p})=\frac{\mathrm{i}e^{\mathrm{i}E\tau_{p}}}{E-\hat{H}_{\ell}}\phi_{\ell m}(r,\tau_{p})\bigr{|}_{R_{\mathrm{I}}},$$

$${\hat{H}}_{\ell }=-\frac{{\partial }_r^2}{2}+\frac{\ell (\ell +1)}{2r^2}+U(r)-\mathrm{i}{V}_{\text{im}}(r).$$

$$({\hat{H}}_{\ell }^{\prime }-\mathrm{i}{\partial }_t)G_{\ell }(r,t)=\delta (t)\delta (r-r^{\prime }).$$

$$\delta a_{\text{I},\ell m}(k,{\tau }_p,T)$$

$$h(t)={\sin }^2(\pi t/T)$$

$$A_x(t)=\frac{E_0}{\omega }{\sin }^2\left(\frac{\omega t}{2n_c}\right)\sin (\omega t),$$

$$A_y(t)=\frac{E_0}{\omega }{\sin }^2\left(\frac{\omega t}{2n_c}\right)\cos (\omega t).$$

$$A_z(t)=\frac{E_0}{\omega }{\sin }^2\left(\frac{\omega t}{2n_c}\right)\sin (\omega t)$$

$$A_x(t)$$

$$A_y(t)$$

$$⟨\ddot{z}(t)⟩=-⟨\Psi (t)\left|\frac{\partial U(r)}{\partial z}\right|\Psi (t)⟩-E_z(t),$$

$$⟨{\ddot{d}}_{\Omega }⟩=-{\int }_0^{\infty }W(t)⟨\ddot{z}(t)⟩e^{-\mathrm{i}\Omega t}dt.$$

$${\left|⟨{\ddot{\mathbf{d}}}_{\Omega }⟩\right|}^2={\left|⟨{\ddot{x}}_{\Omega }⟩\right|}^2+{\left|⟨{\ddot{y}}_{\Omega }⟩\right|}^2$$

$$⟨\ddot{x}(t)+\mathrm{i}\ddot{y}(t)⟩=-⟨\Psi (t)\left|\frac{\partial U(r)}{\partial x}+\mathrm{i}\frac{\partial U(r)}{\partial y}\right|\Psi (t)⟩-(E_x(t)+\mathrm{i}{E}_y(t)).$$

Full text

††journal: Computer Physics Communications

NEW VERSION PROGRAM SUMMARY

Manuscript Title: Qprop with faster calculation of photoelectron spectra

Authors: Vasily Tulsky, Dieter Bauer

Program Title: Qprop

Journal Reference:

*Catalogue identifier:

Licensing provisions:* GNU General Public License, version 3

Programming language: C++

Computer: x86_64

Operating system: Linux

RAM: The memory requirement depends on the desired resolution of the photoelectron spectrum: examples provided with the package need about 1 GB of RAM (with the exception long-wavelength that needs 3.5GB of RAM). If no photoelectron spectrum is produced (hhg example) then only a few MB are necessary.

Number of processors used: Calculation of PES supports parallelization on up to $N_{k}$ processes equal to the number of radial momenta of interest.

Keywords: Time-dependent Schrödinger equation, strong-field ionization, photoelectron spectrum, time-dependent surface flux method, two-color fields.

*External routines/libraries: * GNU Scientific Library, Open MPI (optional).

Catalogue identifier of previous version: ADXB_v2_0

Journal reference of previous version: Comput. Phys. Comm. 207(2016) 452-463

Does the new version supersede the previous version?: Fully supports the functionality of Qprop 2.0.

Nature of problem: Efficient calculation of PES for typical strong-field and attosecond physics ionization scenarios.

Solution method: The time-dependent Schrödinger equation is solved by propagating the electronic wavefunction using a Crank-Nicolson propagator. The wavefunction is represented by an expansion in spherical harmonics. The t-SURFF method in combination with i-SURFV is used to calculate PES.

Reasons for the new version: The i-SURFV method is employed to speed up the calculation of PES.

Summary of revisions: The i-SURFV method is implemented. A set of examples is provided.

Restrictions: The atomic potential needs to be of finite range in case of t-SURFF/i-SURFV usage (i.e., the Coulomb tail is truncated at sufficiently large distances). The laser-matter interaction is described in dipole approximation and velocity gauge.

Additional comments: For additional information see www.qprop.de

Running time: Depends on the laser configuration and on the resolution of PES. Most examples require between 7 to 35 minutes. The longest needs about 1.5 hours.

1 Introduction

Intense-laser-matter experiments brought forward many surprising results that were inaccessible to conventional perturbative theoretical approaches (see, e.g., [1]). As a consequence, new but less rigorous or semi-classical methods have been developed [2, 3, 4, 5] that, however, need to be tested against numerical ab initio solutions. Already the solution of the time-dependent Schrödinger equation (TDSE) for a single active electron in an effective atomic potential and in the presence of a classical, strong laser field can be a demanding computational task [6, 7, 8, 9].

The present paper is devoted to a revised version of Qprop — a position-space TDSE-solver for a single active electron bound in a spherically symmetric potential and subject to an external, time-dependent, space-homogeneous electric field (representing the laser field in dipole approximation). Qprop was introduced in Ref. [7]. A revised version of Qprop employing the time-dependent surface flux method (t-SURFF) for the calculation of PES as proposed in Ref. [10] was published in Ref. [9]. Using t-SURFF, the PES are calculated from the probability flux through a surface located sufficiently far away from the effective range of the binding potential. The time interval over which the flux is captured is limited by the simulation time. Hence, those components of the electronic wavefunction that represent the slowest electrons of interest should reach the t-SURFF surface during the simulation time. In practice, that means that the simulation time might be many times the actual laser pulse duration, in particular for the simulation of ultra-short pulse experiments. In this paper, we introduce Qprop 3.0, where this post-pulse propagation just to capture the slow electrons is avoided using the “trick” proposed in [8] called i-SURFV: once the laser field is off, the evolution of the system is described by a time-independent Hamiltonian, and the contribution to the surface flux after the pulse up to infinity can be calculated in a single step. Refining the formulas used in Qprop 2.0 [9], it is possible to reduce this evaluation to an action of a non-local operator.

The paper is organized as follows. Section 2 contains the mathematical formulation of the upgraded version of t-SURFF (i.e., i-SURFV) that is implemented in Qprop 3.0. In Section 3, the most important functions and data structures are described. Section 4 contains examples. Two examples were already in the Qprop 2.0 paper [9], thus demonstrating nicely the improvement in performance using i-SURFV. A few more demo configurations that may serve as useful templates for a user have been added.

Atomic units $\hbar=|e|=m_{e}=4\pi\epsilon_{0}=1$ are used throughout the paper unless other units are explicitly given.

2 Theoretical basis for i-SURFV

2.1 Hamiltonian and wavefunction

We consider a single active electron, initially bound by the atomic potential, under the influence of a laser field. This system is described by the TDSE

[TABLE]

with the Hamiltonian in velocity gauge

[TABLE]

Here, $U(r)$ is the binding potential of the atom, ${\bf A}(t)$ is the vector potential in dipole approximation (i.e., the electric field is ${\bf E}(t)=-\partial_{t}{\bf A}(t)$), and $V_{\textrm{\scriptsize{ im}}}$ is the imaginary potential which plays the role of an absorber to exclude unphysical reflections in the wavefunction off the numerical boundary. The purely time-dependent term $\sim{\bf A}^{2}(t)$ (that arises from the minimum coupling term $(\hat{{\bf p}}+{\bf A}(t))^{2}/2$) has been transformed away. Due to the spherical symmetry of the potential $U(r)$, it is convenient to expand the wavefunction in spherical harmonics,

[TABLE]

The upper limit for the orbital angular momentum quantum number $\ell$ is finite in numerical calculations, i.e.,

[TABLE]

$N_{\ell}$ is defined by the user. If the polarization of the laser is chosen linear along the $z$ axis and the magnetic quantum number $m_{0}$ of the initial state is well defined, $m_{0}$ is conserved during time propagation so that only $m_{0}$ contributes to the sum over $m$,

[TABLE]

We discretize time and the radial coordinate in units of $\Delta t$ and $\Delta r$, respectively. The radial grid is of size $R_{\mathrm{max}}=R_{\mathrm{I}}+\alpha+W_{\mathrm{im}}$ where $R_{\mathrm{I}}$ is the position of the flux-capturing surface for t-SURFF, $\alpha=E_{\mathrm{max}}/\omega^{2}$ is the quiver amplitude in the chosen laser field of electric field maximum $E_{\mathrm{max}}$, and $W_{\mathrm{im}}$ is the width of the imaginary potential that has the form

[TABLE]

with $R_{\mathrm{im}}=R_{\mathrm{I}}+\alpha$, and $V_{\mathrm{im},\mathrm{max}}=100$ by default. The spherically symmetric potential $U(r)$ can be defined by the user. In the following, we use a potential of the residual atom that is hydrogenic but is switched to linear at $r>R_{\mathrm{co}}$ and is off after reaching zero at $r=2R_{\mathrm{co}}$:

[TABLE]

The cutoff radius $R_{\mathrm{co}}$ is defined by the user in the initial.param file. The t-SURFF/i-SURFV method for the calculation of PES requires that $R_{\mathrm{I}}>2R_{\mathrm{co}}$. The shapes of the binding and imaginary potentials are defined in potentials.hh.

2.2 t-SURFF and i-SURFV

Before we move on to the i-SURFV method, let us briefly review the basics of t-SURFF [10] (see [9] for t-SURFF in the context of Qprop). The PES amplitudes $a_{\textmd{\tiny{I}}}({\bf k},T)$ at time $T$ after the laser pulse are approximated in t-SURFF by projecting the part of the wavefunction that is farther away from the origin than $R_{\mathrm{I}}$ onto Volkov states of the momentum of interest ${\bf k}$,

[TABLE]

Volkov states are plane-wave states in the presence of a laser field [11, 4]. $R_{\mathrm{I}}$ is chosen such that the potential and, hence, the bound states in the region beyond are irrelevant for the effect studied. The time $T$ should be large enough such that an electron with the smallest momentum of interest $k_{\mathrm{min}}$ arrives at the t-SURFF boundary $R_{\mathrm{I}}$, i.e., $T\geq R_{\mathrm{I}}/k_{\mathrm{min}}$. The amplitude (8) can be rewritten in the form of a time integral

[TABLE]

The second term vanishes, as with a properly chosen $R_{\mathrm{I}}$ the initial wavefunction is negligible for $r>R_{\mathrm{I}}$. Using the TDSE (1) for $|\Psi(t)\rangle$ and $-\mathrm{i}\partial_{t}\langle{\bf k}(t)|=\hat{H}_{0}\langle{\bf k}(t)|$ for $\langle{\bf k}(t)|$, where $\hat{H}_{0}=-\Delta/2-\mathrm{i}{\bf A}(t)\cdot\nabla$ coincides with $\hat{H}$ in the region $r>R_{\mathrm{I}}$, the remaining term can be written as

[TABLE]

At this point, let us split the integral over time into parts before and after the end of the laser pulse $\tau_{p}$,

[TABLE]

The first term $a_{\textmd{\tiny{I}}}({\bf k},\tau_{p})$ is treated in Qprop numerically as it was described in [9] while the second, field-free term $\delta a_{\textmd{\tiny{I}}}({\bf k},\tau_{p},T)$ can be calculated without expensive post-propagation of the wavefunction after the laser pulse. This avoidance of post-propagation up to the time where the slowest electrons of interest passed the t-SURFF boundary is the essence of i-SURFV as proposed in Ref. [8] and the core advantage of the new version of Qprop 3.0 over version 2.0.

Let us rewrite the field-free part in the matrix element (11) as

[TABLE]

with $\psi^{*}_{\bf k}({\bf r},t)=\langle{\bf k}(t)|{\bf r}\rangle$ and $\Psi({\bf r},t)=\langle{\bf r}|\Psi(t)\rangle$. Integrating the last term by parts, we obtain

[TABLE]

The (complex conjugated) plane-wave final-momentum state in position space at times $t\geq\tau_{p}$ can be expressed as

[TABLE]

where $j_{\ell}(kr)$ are spherical Bessel functions, and $\Psi({\bf r},t)$ is represented as (4) or (5). Expanding

[TABLE]

one finds

[TABLE]

with $kj^{\prime}_{\ell}(kr)=\partial_{r}j_{\ell}(kr)$. The formula $xj^{\prime}_{\ell}(x)=-xj_{\ell+1}(x)+\ell j_{\ell}(x)$ is used in the program code to calculate the derivatives of the spherical Bessel functions.

The time integral can be evaluated in the limit $T\rightarrow\infty$ using the explicit form of the time evolution operator,

[TABLE]

where the time-independent Hamiltonian has the form

[TABLE]

The operator $\hat{G}_{\ell}(E)=(\hat{H}_{\ell}-E)^{-1}$ is the Green’s function for the radial Schrödinger equation, i.e., the energy representation of the solution to

[TABLE]

The related function in the i-SURFV code is thus named gfunc. The contribution of the upper limit in (17) vanishes due to the presence of an absorbing, imaginary potential. The final expression for the field-free part of the full amplitude can be thus written as

[TABLE]

where $E_{k}=k^{2}/2$. The integral over an infinite time has been converted to a single application of the operator $\left[E_{k}-\hat{H}_{\ell}\right]^{-1}$ (per $k$ and $\ell$), which is the essence of the infinite-time version of t-SURFF, i.e., i-SURFV. How the application of such an inverse operator is actually implemented needs not to be detailed here, as it is similar to the application of the spectral window operator explained in detail in [7]. The applicability of the i-SURFV “trick” is discussed in section 4 by several example configurations included in the Qprop 3.0 distribution.

3 Qprop** general structure**

Parameters and flags

In the current version, all parameters defining coordinate, momentum, and time grids, the potentials, and the laser are moved to the *param files. Most of the flags that allow to switch between different methods or to turn on and off the generation or the storage of specific output are also put into those files. Thus, it is no longer necessary to touch *.cc files for a wide range of problems. All parameters and flags are commented so that their function should become very clear while going through the examples in section 4.

Functions and classes

The core parts of the current Qprop version are as follows.

The real or imaginary time propagation by a single timestep $\Psi(t)\rightarrow\Psi(t+\Delta t)$ is performed in the member function propagate of class wavefunction, which is described in the first Qprop paper [7].

• 2.

The class tsurffSaveWF was designed for saving $\phi_{\ell m}(R_{\mathrm{I}},t)$ and $\partial_{r}\phi_{\ell m}(r,t)|_{r=R_{\mathrm{I}}}$ required for t-SURFF. For i-SURFV, an additional class tsurffSave_full_WF saving $\phi_{\ell m}(r,t)$ was added. The data files end with *.raw.

• 3.

The Green’s function (17) is applied in gfunc, calculating $(\hat{H}_{\ell}-E)^{-1}\phi_{\ell m}(r,t)|_{r=R_{\mathrm{I}}}$ and $\partial_{r}(\hat{H}_{\ell}-E)^{-1}\phi_{\ell m}(r,t)|_{r=R_{\mathrm{I}}}$ from $\phi_{\ell m}(r,t)$ for given $E$.

• 4.

The calculation of a PES using t-SURFF and i-SURFV is performed with the help of the class tsurffSpectrum.

Surface flux output format

In the Qprop 2.0 paper [9], two possible expansions of the amplitudes $a_{\textmd{\tiny{I}}}({\bf k})$ were introduced. If the expansion in angles (according eq. (21) in [9]) is chosen, the data generated by tsurffSpectrum ((29) in [9]) in the output files tsurff_polar$i_{p}$.dat111Here, $i_{p}$ refers to the index of the process that generated the file ($i_{p}=0$ if the non-parallelized version is used). is formatted as

[TABLE]

.

For a linearly polarized laser pulse along the $z$ axis, the spectrum does not depend on the azimuthal angle $\varphi$, and the calculation is performed for $\varphi=0$ only. For a laser pulse polarized in the $xy$ plane, the user might be solely interested in the spectrum in the polarization plane. For that purpose, only data for $\theta=\pi/2$ is generated if $N_{\theta}=1$ is chosen.

If a “complete” expansion in spherical harmonics ((30) in [9]) is desired, data with partial amplitudes according eq. (35) in [9] are stored in tsurff_partial$i_{p}$.dat files in the format

[TABLE]

where $\tilde{\ell}=\ell$ with $N_{\tilde{\ell}}=N_{\ell}-1$ in the case of expansion (5), and $\tilde{\ell}=(\ell+1)\ell+m$ with $N_{\tilde{\ell}}=N_{\ell}^{2}-1$ in the general case (4).

One should note that the complete expansion involves the calculation of Clebsch-Gordan coefficients in terms of Wigner 3$j$ symbols, which are evaluated using the GNU Scientific Library (GSL). However, the respective GSL routine gsl_sf_coupling_3j appears to have acceptable precision only for relatively small $\ell\lesssim 50$ ( see, e.g., Fig. 1 in [12]). Nevertheless, for small $\ell$ the complete expansion can be used to identify dominating angular momenta as a function of energy.

Complexity scaling

The main steps that are required for the calculation of PES with Qprop scale as follows with respect to the number of timesteps, radial grid points, orbital momenta, photoelectron energies and angles:

Imaginary-time propagation: $N_{R}N_{t}$

• 2.

Real-time propagation: $N_{\tilde{\ell}}N_{R}N_{t}$

• 3.

t-SURFF or i-SURFV: $N_{\tilde{\ell}}N_{t}N_{E}N_{\theta}N_{\varphi}$

• 4.

Green’s function “trick” for i-SURFV: $N_{\tilde{\ell}}N_{R}N_{E}$

4 Examples

The quickest way to run an example is to go to its folder and launch the do_all bash script by typing ./do_all.sh in the terminal. Alternatively, one may make and execute the programs for imaginary time propagation (imag_prop), real-time propagation (real_prop) and t-SURFF (tsurff or tsurff_mpi for an MPI-parallelized version) manually. These programs have been slightly revised and renamed in Qprop 3.0. Table 1 shows the old program names in Qprop 2.0, the new names in Qprop 3.0, together with the program task. The user may edit the tsurff.param file and set the variable tsurff-version equal to 1 for t-SURFF or equal to 2 for i-SURFV. If i-SURFV is chosen, it is required to make and execute isurfv after the real-time propagation (this step is included in ./do_all.sh). This generates the data used to calculate the laser-free contribution to the PES according to eq. (2.2). Already in Qprop 2.0, half a Hanning window

[TABLE]

was by default multiplied to the integrand for the $T/2<t<T$ interval of the time integral (13). We keep this feature in Qprop 3.0 if t-SURFF is used.

All parameters defining the $r\ell$-grid, the laser pulse and the momentum grid for the PES are defined in the *.param files. Their precise meaning can be found in the comments there.

All plots shown in this paper were produced using Python scripts that are provided with Qprop 3.0. A brief guide to them can be found at the end of the paper in section 5.

4.1 Speeding up a previous Qprop 2.0 example

Time222For this and other examples, the estimated run times are calculated according [imag_prop + real_prop + isurfv time] + ([tsurff_mpi average time] $\times$ 4 processes used) if the i-SURFV method is used. The compilation time is not counted. Calculations are performed on a desktop PC with Intel(R) Core(TM) i5-6500 CPU.: 6 min + (1 min $\times$ 4).

Memory333The maximum RAM space required is determined by the surface flux that needs to be loaded to calculate the PES with t-SURFF. For the default parameters it is 16 bytes $\times$ 4 (or 3 in the linear-polarization case) $\times N_{\tilde{\ell}}N_{k}N_{\theta}N_{\varphi}$.: 229 MB $\times$ 4.

In the previous Qprop version, t-SURFF required lengthy post-propagation in real time after the laser pulse is over in order to obtain the yield at low momenta. The time interval necessary to capture the PES down to momentum $k_{\mathrm{min}}$ was estimated as $t_{\mathrm{min}}=R_{\mathrm{I}}/k_{\mathrm{min}}$. In contrast, the i-SURFV real-time propagation up to the end of the pulse plus one additional step covering $t\to\infty$ is independent of the lowest momentum of interest. The difference between the total calculation times for the “old” t-SURFF and the new i-SURFV is especially pronounced for short pulses, which are, in fact, often used in modern experiments. To demonstrate this, we revisit an example from Qprop 2.0 [9]: a hydrogen atom in a very short ($n_{c}=2$ cycles) circularly polarized pulse444Note that in the new version of Qprop linear or circular laser polarizations are chosen in initial.param where the value of qprop-dim is set to $34$ for linear polarization (along $z$ axis) or $44$ for circular polarization (in $xy$ plane). No .cc source files need to be modified anymore. with $\omega=0.114$ ($\lambda=400$ nm), $E_{0}=0.0534$ ($\textrm{I}=2\cdot 10^{14}$ Wcm-2*) described by the vector potential ${\bf A}(t)=(A_{x}(t),A_{y}(t),0)$ (the shape of which is defined in potentials.hh),

[TABLE]

The corresponding folder is named attoclock, eponymous to the experiments with such pulses [13]. Compare Fig. 1(a) with Fig. 1(b) that represent the results obtained with t-SURFF and i-SURFV, respectively. The real-time propagation for i-SURFV is equal to the pulse duration and takes 2205 timesteps. It is significantly smaller than the total time required for t-SURFF with $k_{\mathrm{min}}=0.1$, which takes 22205 timesteps. Fig. 1(c) shows the spectra cut along $250^{\circ}$ (and $70^{\circ}$, represented by negative $k$ values). A discrepancy on the logarithmic scale is visible in the high-energy region close to the noise level. This could be improved by increasing the values for $R_{\mathrm{I}}$ and $W_{\mathrm{im}}$, i.e., R-tsurff and imag-width in propagate.param. The low-energy part of the i-SURFV spectrum has a ring-shaped, sharp peak. This artifact is not caused by i-SURFV and would appear in t-SURFF as well if a longer laser-free post-propagation were performed (i.e., if a smaller $k_{\mathrm{min}}$ were chosen). The origin of this structure is the truncation of the potential, and it can be shifted to even lower energies by increasing the cutoff radius $R_{\mathrm{co}}$ (i.e., pot-cutoff in initial.param) or by choosing a smoother potential shape in the intermediate domain $R_{\mathrm{co}}<r<2R_{\mathrm{co}}$. As noted above, the potential should vanish at the flux-capturing surface and beyond so that Volkov states can be used in (8).

Besides the artificial sharp peak, an unphysical oscillatory behavior is present in the low-energy region of Fig. 1(b). This was also observed in [8] and, in fact, is not due to the i-SURFV trick itself but because of the imaginary potential. Upon increasing $W_{\mathrm{im}}$ for the final i-SURFV step according (17), those rings go to lower energies and decrease in magnitude. For that reason a parameter isurfv-imag-width-factor is added to the tsurff.param file.

Figure 2 illustrates the disappearance of these oscillations in the case of linear polarization with the pulse defined according (22). Further, $k_{\mathrm{min}}=0.001$, $N_{k}=2000$, $N_{\theta}=101$, $N_{\varphi}=1$ were set in tsurff.param. Note that only the sharp peak around $k=0.056$ caused by the truncation of the potential survives, and the oscillations due to the imaginary potential are shifted to 10 times smaller momentum.

Alternative ways to deal with the absorber-produced artifacts exist. For instance, a special shape of $V_{\mathrm{im}}(r)$ to supress reflections and transitions might be chosen (see, e.g., [14]) and defined accordingly in potentials.hh. However, this approach is energy-dependent. Another way is to employ complex scaling in the absorption region instead of an imaginary potential [15, 16]. Complex scaling is not yet implemented in Qprop though.

Long wavelength

Time: 46 min + (54 min $\times$ 4).

Memory: 879 MB $\times$ 4.

We now reconsider the most resource-greedy example in Ref. [9]: the calculation of the PES for atomic hydrogen subject to a 6-cycles laser pulse

[TABLE]

with long wavelength $\lambda=2000$ nm and intensity $I=10^{14}$ Wcm*-2*. The sources are located in the folder long-wavelength. In Ref. [9], the focus was on the $10U_{\mathrm{p}}$ rescattering regime ($U_{\mathrm{p}}=E_{0}^{2}/4\omega^{2}$ is the ponderomotive energy). The low-energy part of the PES was not of interest so that post-propagation after the laser pulse was not an issue. Now we apply i-SURFV so that we can examine the low-energy part for this setup, completing the description of the entire spectrum. To speed calculations up, the upper limit $N_{\ell}$ for the angular momenta is decreased to $150$, which is possible because we do not need to describe the high-energy electrons with high accuracy anymore. In fact, $N_{\ell}=150>(I_{\mathrm{p}}+2U_{\mathrm{p}})/\omega$ was chosen instead of $(I_{\mathrm{p}}+10U_{\mathrm{p}})/\omega=623$ required in Ref. [9]. Figure 3 shows spectra obtained with t-SURFF and i-SURFV. Here, $k_{\mathrm{min}}$ for t-SURFF was chosen 0.1 instead of 0.5 in [9]. This slows down the “old” t-SURFF approach because of the longer post-propagation time (93070 timesteps in total). The photoelectron yield near zero momentum is still not converged. Instead, i-SURFV requires just 33070 timesteps and is capable of capturing the low-energy yield.

4.2 Arbitrary two-color pulses

Two-color laser pulses are frequently used in experiments because they allow for additional control of the electron dynamics. For instance, the net current in a plasma after the interaction with a two-color laser pulse can be optimized, which is important for the plasma-induced generation of terahertz radiation. Two-color setups were simulated already with previous versions of Qprop in Refs. [17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. Here, we include illustrative examples that can be easily customized. The corresponding folders are named vortex and lissajous555As a test, we additionally reproduced the single-color example attoclock within the two-color example folder and noticed no difference in the speed of the calculations.. The laser pulse with vector potential ${\bf A}(t)=(A_{x}(t),A_{y}(t),0)$ now reads

[TABLE]

The pulse parameters are specified in the propagate.param file. The matching of the parameter names to the variables appearing in (4.2) is shown in Table 2. Although common wavelengths in the majority of modern two-color experiments are about 800 nm, we take for convenience $\lambda\lesssim 400$ nm in the examples here in order to keep the calculation times short.

Vortex

Time: 25 min + (7 min $\times$ 4).

RAM: 206 MB $\times$ 4.

In the example located in the vortex folder, the laser field consists of two sequential, counter-rotating circularly polarized pulses with $\omega_{1}:\omega_{2}=3:5$ ($\lambda_{1}=400$ nm, $\lambda_{2}=240$ nm), equal intensities $I=5\cdot 10^{13}$ Wcm*-2* and equal total time duration $\tau=8$ fs. Note that all field component amplitudes $E_{ij}$ have to be non-negative, and the opposite direction of the $\omega_{2}$ field is achieved by adding $\pi$ to the phase $\phi_{2y}$ in the $y$ component. For simplicity, the hydrogen atom with the ionization potential $I_{\mathrm{p}}=0.5$ is used as a target. The ionization potential may be overcome by five $\omega_{1}=0.114$ photons or by three $\omega_{2}=0.190$ photons. Thus, photoelectrons with energy $k^{2}/2$ are represented by a superposition of wavefunctions with $Y_{55}(\theta,\varphi)\propto e^{5\mathrm{i}\varphi}$ and $Y_{3-3}(\theta,\varphi)\propto e^{-3\mathrm{i}\varphi}$ angular dependencies. Additionally, the delay $\tau$ between the pulses creates a phase factor $e^{-\mathrm{i}k^{2}\tau/2}$ so that the photoelectron yield behaves as $k|a({\bf k})|^{2}\propto A+B\cos(8\varphi-k^{2}\tau/2+\alpha)$, resulting in an eight-fold vortex-like structure (see Fig. 4(a)) with a torsion defined by $\tau$ (see the tilt in Fig. 4(b)). Due to the finite duration of the pulses, frequencies with ratios other than $3:5$ may contribute so that the total number of nodes differs for energies away from the maximum yield. Such kind of spectra were observed in Refs. [27, 28, 29, 30] in which one may also find more detailed explanations of the underlying physics.

Lissajous

Time: 14 min + (4 min $\times$ 4).

RAM: 206 MB $\times$ 4.

The example in the lissajous folder refers to another type of laser pulses that has been used in recent theoretical and experimental studies, namely two overlapping, orthogonally linearly polarized pulses with frequencies of ratio $n:m$ as, e.g., in Refs. [17, 18, 24]. In our example, we take $\lambda_{1}=400$ nm and $\lambda_{2}=200$ nm fields of equal intensity $I=10^{14}$ Wcm*-2* and equal time duration $10.7$ fs. The laser field vectors ${\bf E}(t)$ (or ${\bf A}(t)$) form Lissajous curves in the polarization plane (see Fig. 5(b,c)). Such pulse geometries allow to separate effects due to scattering off the parent ion from intracycle interferences.

Two-color, colinear polarization

Two both along the $z$-axis colinearly polarized laser pulses can be simulated using qprop-dim $=34$, which allows to use the wavefunction expansion in spherical harmonics with fixed magnetic quantum number $m_{0}$ (5). In that case the $x$ components are automatically taken as $z$ components in the code.

Such two-color, colineraly polarized pulses are, for instance, common in experiments concerned with plasma-induced terahertz radiation [31] where they act as ionizing fields that prepare the initial plasma [32]. Another example is the addition of a weak $2\omega$ field to study the dependence of PES on the delay between strong $\omega$ and weak $2\omega$ field. Besides the additional control knob to steer the electron dynamics, such “phase-of-the-phase” [33, 34, 35, 22, 36, 37, 38, 20, 39] studies allow to reveal the coherently produced contributions to the PES that may otherwise be buried under incoherent, e.g., thermal or scattered, electrons.

4.3 High-order harmonic generation

Time: 7 min.

Besides PES, strong-field physicists are often interested in the radiation emitted by laser-illuminated targets. It is well known that multiples of the incident laser frequency are emitted [40, 41] due to so-called high-harmonic generation (HHG). Computationally, the calculation of high-harmonic spectra is much simpler than the simulation of PES, at least as long as a semi-classical description and the single-atom response are sufficient. In that case only the dipole acceleration needs to be calculated, and the PES-related t-SURFF part can be omitted.

We include an example for the calculation of a HHG spectrum in the folder hhg. There, a linearly polarized pulse with an $n_{c}-4=6$-cycle flat-top part and 2-cycle $\sin^{2}$-shape up and down rampings impinges on a hydrogen atom. Wavelength and intensity are $\lambda=800$ nm and $I=10^{14}$ Wcm*-2*, respectively. Figure 6 shows the expectation value of the acceleration as a function of time, calculated as

[TABLE]

and the HHG spectrum $\left|\langle\ddot{d}_{\Omega}\rangle\right|^{2}$ with

[TABLE]

The window function $W(t)$ was chosen $\sin^{2}$ in the flat-top region of the laser pulse and zero outside. That allows to obtain very pronounced HHG peaks with a high signal-to-noise ratio. The HHG spectrum consists of sharp peaks at odd multiples of the incident laser frequency, forming the so-called plateau up to the cutoff $\Omega=I_{\mathrm{p}}+3.17U_{\mathrm{p}}$ [42, 41] and a fast decrease thereafter.

For a laser pulse polarized in the $xy$ plane, the HHG spectrum is calculated according

[TABLE]

where first and second term are calculated separately from the Fourier-transforms of the real and imaginary parts of

[TABLE]

The ellipticity of a certain harmonic can be obtained from the phase difference of the Fourier transforms at the respective $\Omega$ value.

The reader may, for instance, reconsider the lissajous example, switch the generate-hhg-data trigger from 0 to 1, the generate-tsurff-data trigger from 1 to 0 and set pot-cutoff larger than the grid size to avoid a potential truncation (required for t-SURFF but not for HHG). In such an orthogonal, two-color $\omega$-$2\omega$ setup the dipole acceleration expectation value in $y$ direction consists of even harmonics of the fundamental frequency $\omega$ (in $x$ direction). Therefore, the sum (26) yields a HHG spectrum with both odd and even multiples of $\omega$. For a more detailed study and applications see, e.g., Refs. [43, 44, 45, 46, 47].

5 Plotting guide

For the user’s convenience, scripts written in Python that were used to visualize the Qprop 3.0 generated data are added to the package. They are located in scr/plots.

Photoelectron distributions can be plotted using plot_pes.py. Leave the desired filename uncommented in the upper section of it, choose the number of angles and the polarization as =’xz’ for linear or =’xy’ otherwise.

Figure 1(a) was produced with polar_canvas=1, Fig. 1(c) with plot_type=’1D’, Fig. 2(a) was based on the same script with polar_canvas=0 and cartesian=1. Figure 4(b) shows a spectrum plotted with respect to energy, hence wrt_energy=1, and a canvas format polar_canvas=0, cartesian=0. Laser field plots in Fig. 5(b,c) are generated with the script plot_laser.py. Plots in Fig. 6 are created with a slight modification of the script plot_hhg.py.

Acknowledgment

This work was supported by the project BA 2190/10 of the German Science Foundation (DFG).