Fiend -- Finite Element Quantum Dynamics

TL;DR

Fiend is a Python-based simulation package for solving the 3D time-dependent Schrödinger equation in cylindrically symmetric systems, enabling detailed studies of electron dynamics under complex electromagnetic fields with advanced numerical methods.

Contribution

It introduces a finite element method-based solver with full minimal coupling for inhomogeneous vector potentials, supporting flexible and accurate quantum dynamics simulations.

Findings

Supports inhomogeneous vector potentials beyond dipole approximation

Uses FEM for spatial discretization and advanced time-stepping methods

Provides user-friendly API with extensive documentation

Abstract

We present Fiend - a simulation package for three-dimensional single-particle time-dependent Schr\"odinger equation for cylindrically symmetric systems. Fiend has been designed for the simulation of electron dynamics under inhomogeneus vector potentials such as in nanostructures, but it can also be used to study, e.g., nonlinear light-matter interaction in atoms and linear molecules. The light-matter interaction can be included via the minimal coupling principle in its full rigour, beyond the conventional dipole approximation. The underlying spatial discretization is based on the finite element method (FEM), and time-stepping is provided either via the generalized-{\alpha} or Crank-Nicolson methods. The software is written in Python 3.6, and it utilizes state-of-the-art linear algebra and FEM backends for performance-critical tasks. Fiend comes along with an extensive API documentation, a user guide, simulation examples, and allows for easy installation via Docker or the Python Package Index.

Figures

Equations

$$\begin{array}{rl}\hat{H}& =\widehat{T_{\rho }}+\frac{{\hat{p}}_z^2}{2}+V(\hat{\rho },\hat{z})+W(\hat{\rho },\hat{z},{\hat{p}}_{\rho },{\hat{p}}_z,t),\end{array}$$

$$\begin{array}{rl}{\hat{T}}_{\rho }& =-\frac{1}{2}\left(\frac{{\partial }^2}{{\partial }_{\rho }^2}+\frac{1}{2}\frac{\partial }{{\partial }_{\rho }}\right)\\ & \ne \frac{1}{2}{\left[-\mathrm{i}\left(\frac{\partial }{{\partial }_{\rho }}+\frac{1}{2\rho }\right)\right]}^2=\frac{{\hat{p}}_{\rho }^2}{2}.\end{array}$$

$$W_{\mathrm{LG}}=-\frac{\partial f(t)}{\partial t}\hat{z},$$

$$W_{\mathrm{VG}}=f(t){\hat{p}}_z,$$

$$\begin{array}{rl}{W}_{\mathrm{inhomogeneous}}& =\frac{1}{2}f(t){\mathbf{A}}_s\cdot \hat{\mathbf{p}}+\frac{1}{2}f(t)\hat{\mathbf{p}}\cdot {\mathbf{A}}_s\\ & +\frac{1}{2}f(t{)}^2\Vert {\mathbf{A}}_s(\hat{\rho },\hat{z}){\Vert }^2.\end{array}$$

$$-\frac{1}{2\rho }\frac{\partial }{\partial \rho }\left(\rho \frac{\partial {\psi }_k}{\partial \rho }\right)-\frac{1}{2}\frac{{\partial }^2{\psi }_k}{\partial z^2}+V(\rho ,z){\psi }_k=E_k{\psi }_k,$$

$$\underset{\rho \to 0+}{\lim }\left(\rho \frac{\partial \psi }{\partial \rho }\right)=0.$$

$$\begin{array}{rl}& -\frac{1}{2}\underset{{\Omega }_{\text{TI}}}{\int }\sum _{\alpha =\rho ,z}\frac{\partial \chi }{\partial \alpha }\frac{\partial {\psi }_k}{\partial \alpha }\rho \mathrm{d}\rho \mathrm{d}z\\ & +\underset{{\Omega }_{\text{TI}}}{\int }\chi V(\rho ,z){\psi }_k\rho \mathrm{d}\rho \mathrm{d}z\\ & =E_k\underset{{\Omega }_{\text{TI}}}{\int }\chi {\psi }_k\rho \mathrm{d}\rho \mathrm{d}z\forall \chi \in F,\end{array}$$

$$\text{innerproduct}\chi \psi =\underset{{\Omega }_{\text{TI}}}{\int }\chi (\rho ,z{)}^{\ast }\psi (\rho ,z)\rho \mathrm{d}\rho \mathrm{d}z.$$

$$F_h=\mathrm{span}{\left\{{\varphi }_i\right\}}_{i=0}^{N-1}\subset F$$

$$(\mathbf{T}+\mathbf{V}){\mathbf{\psi }}_k=\mathbf{S}{E}_k{\mathbf{\psi }}_k.$$

$$S_{ij}$$

$$T_{ij}$$

$$V_{ij}$$

$$\mathrm{i}\mathbf{S}\dot{\mathbf{\psi }}(t)=\left(\mathbf{T}+\mathbf{V}+\mathbf{W}\right)\mathbf{\psi }(t),$$

$$W_{ij}=\text{matrixelement}{\varphi }_i{W}_s\left(\rho ,z,\frac{\partial }{\partial \rho },\frac{\partial }{\partial z},t\right){\varphi }_j$$

$$\begin{array}{rl}\psi (t+\Delta t)=& {\left[\mathbf{S}+\frac{\mathrm{i}\Delta t}{2}\mathbf{H}\left(t+\frac{\Delta t}{2}\right)\right]}^{-1}\\ & \left[\mathbf{S}-\frac{\mathrm{i}\Delta t}{2}\mathbf{H}\left(t+\frac{\Delta t}{2}\right)\right]\psi (t)+\mathcal{O}(\Delta t^2),\end{array}$$

$${\hat{p}}_{\rho }\to -\mathrm{i}\left[{\partial }_{\rho }+\frac{1}{2\rho }-\frac{1}{2}\delta (\mathbf{r}-{\Gamma }_N)n_{\rho }(\mathbf{r})\right],$$

$${\hat{p}}_z\to -\mathrm{i}\left[{\partial }_z-\frac{1}{2}\delta (\mathbf{r}-{\Gamma }_N)n_z(\mathbf{r})\right],\text{ and}$$

$$\frac{1}{2}(\mathbf{A}\cdot \hat{\mathbf{p}}+\hat{\mathbf{p}}\cdot \mathbf{A})\to \frac{1}{2}\left(\mathbf{A}\cdot \hat{\mathbf{p}}+\hat{\mathbf{p}}\cdot \mathbf{A}\right)$$

$$-\frac{\mathrm{i}}{2}\mathbf{A}\cdot \mathbf{n}(\mathbf{r})\delta (\mathbf{r}-{\Gamma }_N),$$

$$\gamma \sum _i{\left|\underset{{\Gamma }_N}{\int }{\chi }_i^{\ast }(\rho ,z)\nabla \psi (\rho ,z,t)\cdot \mathrm{d}\mathbf{S}\right|}^2,\gamma \in \mathbb{R}$$

$$\frac{{\mathbf{\psi }}_k^{†}{\mathbf{H}}_0{\mathbf{\psi }}_k}{{\mathbf{\psi }}_k^{†}\mathbf{S}{\mathbf{\psi }}_k}$$

$$\begin{array}{r}\max \text{CR}(\mathbf{r})={\text{CR}}_{\text{asymp}}[1-\left(1-\frac{{\text{CR}}_{\text{ref}}}{{\text{CR}}_{\text{asymp}}}\right)\\ \times \frac{R_{\text{ref}}}{r}\tanh \left(\frac{r}{R_{\text{ref}}}\right)\exp \left(-\frac{r^2}{R_{\text{trans}}^2}\right)].\end{array}$$

$$\psi (t=0)=\frac{1}{\sqrt{2}}\left({\psi }_{1s}+{\psi }_{2s}\right).$$

$$\left|1-|⟨{\psi }_{\text{exact}}||{\psi }_{\text{Fiend}}⟩{|}^2\right|.$$

$$\left|1-\Vert {\psi }_{\text{Fiend}}{\Vert }^2\right|.$$

$$S(\omega )\propto {\left|\text{expectationvalue}{\ddot{D}}_z(\omega )\right|}^2,$$

$$\text{expectationvalue}{\ddot{D}}_z(t)=\text{matrixelement}\psi (t)-\frac{\partial V}{\partial z}\psi (t).$$

$$\mathbf{A}(\mathbf{r},t)={\mathbf{A}}_s(\mathbf{r})f(t),$$

$${\mathbf{A}}_s(\mathbf{r})=\mathbf{E}(\mathbf{r})=-\nabla U(\mathbf{r}).$$

$$-\nabla \cdot \left[ϵ(\mathbf{r})\nabla U(\mathbf{r})\right]=0,$$

$${ϵ}_{\mathrm{Au}}\nabla u{|}_{\mathbf{r}\to {\Gamma }_{\text{tip-vac}}}-{ϵ}_{\mathrm{vac}}\nabla u{|}_{\mathbf{r}\to {\Gamma }_{\text{vac-tip}}}=0$$

Full text

††journal: Computer Physics Communications

**PROGRAM SUMMARY

** Program Title: Fiend 

Journal Reference:

Catalogue identifier:

Licensing provisions: MIT License

Programming language: Python 3.6

Computer: Tested on x86_64 architecture.

Operating system: Tested on Linux and macOS.

RAM: Simulation dependent: from megabytes to hundreds of gigabytes.

Parallelization: MPI-based parallelization. Thread-based parallelization via BLAS and LAPACK backends.

Classification: 2.2 Spectra, 2.5 Photon Interactions, 4.3 Differential equations, 6.5 Software including Parallel Algorithms

External routines/libraries:

PETSc, SLEPc, FEniCS, HDF5

petsc4py, slepc4py, mpi4py, h5py, numpy, scipy, matplotlib, psutil, mypy, progressbar2

*Nature of problem:

*Solution of time independent and time dependent single active electron Schrödinger equations in cylindrically symmetric systems including interactions with spatially inhomogeneous vector potentials.

*Solution method:

*Finite element discretization of the Schrödinger equation. Time evolution via the generalized-$\alpha$ or Crank-Nicolson methods.

*Restrictions:

*Cylindrically symmetric single active electron systems.

*Unusual features:

*Finite element discretization of the equations allowing inhomogeneous spatial dependence of the vector potential (e.g., plasmon-enhanced fields). Integration with the FEniCS finite element suite. Installable from the Python Package Index. Pre-installed Docker images available.

*Additional comments:

*The source code is also available at https://fiend.solanpaa.fi.

Python package available at https://pypi.org/project/fiend/.

Docker images are available at https://hub.docker.com/r/solanpaa/fiend/.

*Running time:

*From minutes to weeks, depending on the simulation.

1 Introduction

Simulation of three-dimensional (3D) single particle quantum mechanics (QM) is still one of the most used computational approaches in the strong field and attosecond communities. In these fields, the nonlinear interaction of the atomic or molecular electron with the driving laser field requires fast and accurate integration of the 3D time dependent Schrödinger equation (TDSE). Recent applications include, e.g., benchmarking of approximate models [1, 2, 3, 4, 5, 6], study of high-order harmonic generation [7, 8, 9], and the study of photoionization [10, 11, 12, 13, 14].

Recently the strong-field and attosecond communities have turned their attention to related phenomena in nanostructures. There the nanostructure geometry and plasmonic effects cause the electromagnetic field to become extremely inhomogeneous [15] which, in turn, causes significant differences to the traditional ultrafast and strong-field phenomena in atoms and molecules. Recent studies include, e.g., nanostructure-enhanced photoionization of gases [16, 17, 18] and electron emission from nanostructures such as tips [19, 20, 21, 22, 23, 24, 25, 26, 27, 28] and rods [29, 30, 31, 32, 33].

There are already multiple software designed for integrating the three-dimensional time-dependent Schrödinger equation (TDSE). Open-source software include, e.g., Qprop [34], Octopus [35, 36], QnDynCUDA [37], and WavePacket [38], but there are also plenty of other options such as SCID-TDSE [39], tRecX [40, 41], and ALTDSE [42]. Qprop, SCID-TDSE, and QnDynCUDA solve the TDSE using a grid-based representation of the radial coordinates and spherical harmonics for the angular dependence; TRecX and WavePacket support multiple basis sets; Octopus relies on a real-space grid; and ALTDSE requires end-user to provide the matrices in an appropriate eigenbasis.

However, efficient and accurate simulation of ultrafast and strong-field phenomena in nanostructures and nanostructure-enhanced gases requires specialized TDSE-solvers. Most importantly, since the laser electric field has strong spatial inhomogeneity, it is imperative to have position-dependent control of simulation accuracy. In addition, it would be beneficial to have the TDSE-solver directly interface with a solver for the electromagnetic problem.

To this end, we have developed Fiend, a QM simulator based on finite element method (FEM). Fiend provides solvers for the time independent and time-dependent single active electron Schrödinger equations in cylindrically symmetric systems111Although the restriction to cylindrically symmetric systems can be lifted relatively easily if needed. for solving time-independent and time-dependent Schrödinger equations in cylindrically symmetric systems. We have designed Fiend to integrate with the open source FEniCS finite element suite [43, 44, 45, 46, 47, 48, 49] allowing for easy description of complicated system geometries. Being FEM-based, Fiend also allows the use of potentials with integrable singularities.

Fiend is written in Python 3.6, but much of the heavy number-crunching is delegated to well-tested and efficient external libraries. The software is modular and easy to extend. We also provide a comprehensive unit- and integration-test suite to ensure reliability. Moreover, we provide straightforward installation methods either via the Python Package Index [50] or as Docker images [51].

This paper is organized as follows. In Sec. 2 we describe the finite element discretization of the Schrödinger equations and the time-propagation schemes. Section 3 focuses on the design and implementation of Fiend. Section 4 provides a few example simulation, and finally, in Sec. 5 we summarize this paper and possible extensions of the presented software suite.

2 Finite element quantum mechanics

2.1 Systems

Fiend is designed for simulating cylindrically symmetric single (active) electron systems interacting with time and (optionally) position dependent vector potentials. The Hamiltonian operator of these systems can be written as

[TABLE]

where $\hat{\rho}$ is the planar radial coordinate operator, $\hat{z}$ the $z$-coordinate operator, $\hat{T}_{\rho}$ the kinetic energy operator with respect to planar radial motion222 Note that $\hat{T}_{\rho}$ is not proportional to the square of the planar radial momentum operator, i.e.,

\begin{split}\hat{T}_{\rho}&=-\frac{1}{2}\left(\frac{\partial^{2}}{\partial_{\rho}^{2}}+\frac{1}{2}\frac{\!\partial}{\partial_{\rho}}\right)\\ &\neq\frac{1}{2}\left[-{\mathrm{i}\mkern 1.0mu}\left(\frac{\!\partial}{\partial_{\rho}}+\frac{1}{2\rho}\right)\right]^{2}=\frac{\hat{p}_{\rho}^{2}}{2}.\end{split}

(2)

This is a peculiarity of quantum mechanics in cylindrical coordinates [52]., $\hat{p}_{z}$ the $z$-component of the momentum operator, $V$ the static potential, and $W$ the light-matter interaction operator.

Fiend readily supports three types of interaction operators $W$. First, for linearly polarized vector potentials $\mathbf{A}=f(t)\bm{e}_{z}$ in the length gauge the interaction operator is given by

[TABLE]

secondly, in the velocity gauge the interaction operator reads

[TABLE]

and finally for more general cylindrically symmetric inhomogeneous vector potentials $\mathbf{A}=\mathbf{A}_{s}(\rho,z)f(t)$ we provide the interaction operator

[TABLE]

Moreover, Fiend supports all interaction operators of the form $\hat{W}=W(\hat{\rho},\hat{z},\hat{p}_{\rho},\hat{p}_{z},t)$.

The above interaction operators conserve the magnetic quantum number $m$, and by default, we simulate only the $m=0$ subspace. Correspondingly, the coordinate space domain is a two-dimensional slice of the cylindrical coordinates, $\Omega_{\infty}=\left\{\rho\geq 0,z\in\mathbb{R}\right\}$, which we truncate for numerical simulations to a finite domain $\Omega$ [see Fig. 1(a)].

2.2 Time independent Schrödinger equation

On the domain $\Omega_{\text{TI}}$ the time-independent Schrödinger equation (TISE) can be written as

[TABLE]

where $E_{k},\psi_{k}$ is the $k$th eigenpair. This equation is accompanied by continuity boundary condition at $\rho=0$,

[TABLE]

Elsewhere on the boundary $\partial\Omega_{\text{TI}}$ we can choose freely between zero Neumann boundary conditions (ZNBCs) and zero Dirichlet boundary conditions (ZDBCs) – depending on the requirements of our model. For example, in Fig. 1(a) where we have imposed ZNBC on the upper boundary of the simulation domain, $\Gamma_{N}$, and ZDBC on the arc $\Gamma_{D}$.

To discretize TISE, we begin by looking for eigenpairs $E_{k}\in\mathbb{R},\psi_{k}\in F$ of the weak form corresponding to Eq. (6) and the boundary conditions described above. The weak form is given by

[TABLE]

where $F=\left\{\psi\in H^{1}(\Omega_{\text{TI}})\,|\,\psi(\Gamma_{D})=0\right\}$ is the standard Sobolev space on $\Omega_{\text{TI}}$ of real-valued $L^{2}$-integrable functions with $L^{2}$ integrable first derivatives. Notice also that $F$ includes ZDBCs, and the natural inner product is

[TABLE]

Furthermore, we must restrict ourselves to a finite dimensional approximation of $F$. First, the simulation domain $\Omega_{\text{TI}}$ is described using an unstructured triangular mesh which supports arbitrary refinement [see, e.g., Fig. 1(b)]. Next, we construct a basis of continuous low-order Lagrange polynomials $\phi_{i}$, each with compact support on the mesh elements. This basis spans a finite dimensional function space

[TABLE]

where the TISE weak form [Eq. (8)] can be written as a finite dimensional generalized Hermitian eigenvalue equation

[TABLE]

Here $\bm{\psi}_{k}$ is a vector of the real-valued expansion coefficients of the $k$th eigenstate, $E_{k}$ the corresponding eigenvalue, and

[TABLE]

are the overlap, kinetic energy, and potential energy matrices, respectively. Note that since the static potential of the system is implemented only via Eq. (14) which involves integration, we can use any static potential for which Eq. (14) is finite, e.g., the Coulomb potential.

2.3 Time dependent Schrödinger equation

In a similar manner as with TISE, TDSE can be discretized as

[TABLE]

where $\bm{\psi}(t)$ is a vector of the complex-valued expansion coefficients $c_{i}(t)$ of the wave function $\psi(\rho,z,t)=\sum c_{i}(t)\phi_{i}(\rho,z)$. $\bm{S}$, $\bm{T}$, and $\bm{V}$ are constructed the same way as for TISE (except we add an imaginary absorbing potential to $\bm{V}$), and the light-matter interaction matrix is given by

[TABLE]

The basis of the discrete function space for TDSE can, in general, be different from the one used for TISE. This is useful in practice since often the stationary states can be computed in a much smaller simulation domain than needed for an accurate description of the TDSE. We can change the function space by three methods: (1) using a larger simulation domain for the TDSE, (2) refining the mesh according to the requirements of the TDSE simulation, and (3) change the degree of our basis functions according to the required accuracy.

For evolving the discretized state, i.e., the expansion coefficients $\bm{\psi}(t)$ according to Eq. (15), Fiend implements two approximations of the time evolution operator: the Crank-Nicolson (CN) method [53] and the generalized-$\alpha$ method [54].

The Crank-Nicolson method can be written as

[TABLE]

and it is a unitary time-reversible transformation [53] as long as the matrix inversion is performed to a sufficient accuracy.

The generalized-$\alpha$ method was originally developed for integrating the equations of motion arising from fluid dynamics [54], but we have successfully applied it to the TDSE integration. To the best of our knowledge, the $\alpha$-method has not been proven to be unitary nor is it exactly time-reversible. Nevertheless, we have achieved accuracies comparable to the CN method – and in the case of extremely inhomogeneous vector potentials, the generalized-$\alpha$ method seems to be significantly more stable than CN.

2.4 Incorporation of boundary conditions

The continuity boundary condition at $\rho=0$ and ZNBCs are automatically included by dropping the boundary integrals when deriving Eqs. (8) and (13). However, the ZDBCs must be included in the system matrices via the following modifications: the rows and columns of $\bm{S}$, $\bm{T}$, $\bm{V}$, and $\bm{W}$ corresponding to the degrees of freedom (DOFs) on the boundary $\Gamma_{D}$ with ZDBCs set to zero. Only for $\bm{S}$ we must set the corresponding diagonal elements to one to ensure invertibility. Note that it’s crucial to zero out not only the rows but also the corresponding columns of the matrices to retain hermiticity.

Embedding ZNBCs via the vanishing boundary terms when deriving Eqs. (8) and (13) causes a practical issue: Some operators such as $\hat{p}_{z}=-{\mathrm{i}\mkern 1.0mu}\partial_{z}$, $\hat{p}_{\rho}=-{\mathrm{i}\mkern 1.0mu}\left(\partial_{\rho}+\frac{1}{2\rho}\right)$, and $\bm{A}\cdot\bm{\hat{p}}+\bm{\hat{p}}\cdot\bm{A}$ become non-Hermitian. Minor adjustments at the boundaries can remedy this, namely,

[TABLE]

where $\bm{n}(\bm{r})=[n_{\rho}(\bm{r}),n_{z}(\bm{r})]^{T}$ is the outwards facing unit normal on the simulation domain boundary at point $\bm{r}$. This trick has been previously used to obtain the correct Hermitian operators in hyperspherical coordinates [52].

Another aspect to consider is the weak enforcement of ZNBCs compared to ZDBCs. If there was significant electron density on the Neumann boundary $\Gamma_{N}$, it will start to oscillate and eventually violate the ZNBC. This can be remedied when operating with matrix inverses as in the second step of the CN-propagator. We can add an extra error term, e.g., of the form

[TABLE]

to the error estimator in our numerical implementation. This will remove – or at least lessen – the issue arising from the weak enforcement of ZNBCs.

3 Implementation

3.1 Overview

We have implemented finite element discretization of TISE and TDSE in the Python software package Fiend following the recipes of Sec. 2. Fiend is written in Python 3.6 and it utilizes reliable libraries commonly pre-installed in high-performance clusters and supercomputers. For sparse linear algebra we use the Portable, Extensible Toolkit for Scientific Computation (PETSc) [55, 56, 57] and the Scalable Library for Eigenvalue Problem Computations (SLEPc) [58, 59]. Meshing and other standard FEM-related parts are built on top of the components of the FEniCS project [43, 44, 45, 46, 47, 48, 49], and filesystem IO is largely based on HDF5 [60] via h5py [61, 62]. In postprocessing and visualization we use also numpy [63], scipy [64], and matplotlib [65]. Fiend is parallelized using MPI via mpi4py [66, 67, 68].

The PETSc-dependency is slightly complicated as PETSc needs to be compiled either with real number support (x)or with complex number support. Fiend, on the other hand, needs real numbers for TISE and post-processing but complex numbers for propagation. Consequently, the user must install both the real and complex versions of PETSc/SLEPc stacks and switch between these for different simulation steps. This factitious requirement for two different installations of PETSc can be revisited in the future upon the completion of DOLFIN-X [69] and FFC-X [70] including the merge of complex number support in FEniCS [71],

To ease installation, we provide a Docker [51] container with Fiend and all its dependencies preinstalled, see https://hub.docker.com/r/solanpaa/fiend/ for more details. In the same spirit, Fiend is also available in the Python Package Index [50] and can be installed with pip [72].

3.2 Library usage and numerical methods

A typical usage of the Fiend suite is demonstrated in Fig. 2. First we solve the TISE to obtain a set of stationary states. This includes meshing the domain for the time independent simulation and assembling the system matrices (12)–(14) with appropriate boundary conditions (see Sec. 2.4).

The TISE eigenproblem, Eq. (11), is solved (by default) with Rayleigh Quotient Conjugate Gradient (RQCG) method combined with classical Gram-Schmidt orthogonalization of the Krylov subspace basis with adaptive iterative refinement for increased numerical stability. RQCG is a variational method which essentially minimizes the Rayleigh quotient

[TABLE]

of a desired number of orthonormal vectors with respect to the bilinear product induced by $\bm{S}$ [59]. Upon convergence, this corresponds to $k$ smallest real eigenpairs of TISE.

Step 2 (see Fig. 2) is to prepare the discrete description of the TDSE. This includes meshing the domain for the time dependent problem, interpolating stationary states to the new mesh, and also assembly of the system matrices  (12)–(14) and the time-independent part(s) of Eq. (16).

The propagation in Step. 3 (Fig. 2) requires the user to load the environment with complex number PETSc. The time-dependent Hamiltonian is constructed from the matrices computed in the previous step, and the discrete TDSE (15) is solved with PETSc’s propagators. We also provide code templates for easy implementation of new propagators.

Finally, we provide a set of example scripts for postprocessing. These scripts include, e.g., visualization of the integrated electron density in time and temporal shape of the laser field, but there are also scripts to compute more complex observables such as the high-harmonic and photoelectron spectra.

All the numerical methods in Fiend depend on efficient sparse linear algebra operations implemented in PETSc. By default, matrix inversions and the solution of linear equations are carried out with Generalized Minimal RESidual method (GMRES) [73], where the convergence criterion is computed with respect to the norm induced by the inner product in Eq. (9). Unfortunately, the PETSc linear algebra backend does not allow us to build the basis of the Krylov subspace with respect to a custom inner product. However, we have found the basis built with the standard Hermitian inner product to be sufficient performance-wise. Direct solvers such as SuperLU_DIST [74, 75, 76] and MUMPS [77, 78] are supported only with matrices that can be explicitly constructed without too high cost.

3.3 Note on meshing

We provide a way to generate meshes with arbitrary refinement in the coordinate space. The user should supply a function returning the maximum allowed cell circumradius333Radius of the smallest circle enclosing the given triangular cell. at coordinate $\bm{r}$, and the mesh can be refined until all cells of the mesh have circumradius below the desired one.

By default, the maximum cell circumradius is given by [79]

[TABLE]

For $r<R_{\text{ref}}$ the cell circumradii are below $\text{CR}_{\text{ref}}$, and as $r\to R_{\text{trans}}$, the maximum cell circumradius increases monotonically to $\text{CR}_{\text{asymp}}$.

4 Examples and test cases

4.1 Field-free propagation

To assess the stability and accuracy of the numerical methods implemented in Fiend, we compute a field-free propagation of an electron prepared in a superposition of Hydrogen 1s and 2s states,

[TABLE]

This state is propagated up to $100\text{\,}\mathrm{a.u.}$ of time. The simulation domain up to $R=$ $30\text{\,}\mathrm{a.u.}$ is meshed using the default refinement introduced in Sec. 3.3 with $R_{\text{ref}}=$ $4\text{\,}\mathrm{a.u.}$, $R_{\text{ref}}=$ $10\text{\,}\mathrm{a.u.}$, $\text{CR}_{\text{ref}}=$ $0.01\text{\,}\mathrm{a.u.}$ and $\text{CR}_{\text{asymp}}=$ $0.5\text{\,}\mathrm{a.u.}$ We assess the accuracy of the numerical solution by comparing the simulated state to the exact result:

[TABLE]

The stability is investigated by evaluating how much the norm of the state differs from unity, i.e.,

[TABLE]

Both the $\alpha$-propagator and the CN propagator reach an accuracy of 0.00524 % for the overlap with the time-step $0.05\text{\,}\mathrm{a.u.}$, and the norm of the simulated state deviates from unity only by $3.22\cdot 10^{-5}\,\%$. This demonstrates the applicability of the propagators for the simulation of FE-discretized TDSE in cylindrical coordinates.

4.2 High-harmonic generation

Next we demonstrate the applicability of Fiend to simulate high-order harmonic generation (HHG). According to the three-step model [80, 81], an atomic electron is excited to continuum and driven to oscillate around the atom [82]. Upon recombination with the parent ion, the excess energy is released as high-energy photons.

Full 3D simulation of HHG via TDSE is a difficult task due to the long extent of the electron wave function when driven by a short laser pulse. In order to obtain a decent description of the single-atom response for HHG, a large simulation domain is needed. We setup a simulation domain of radius $800\text{\,}\mathrm{a.u.}$ with a more refined mesh at the origin and gradually sparser mesh for larger radii.

The intensity spectrum of the single-atom response is proportional to the square of the Fourier transform of the dipole acceleration [83],

[TABLE]

where we compute the dipole acceleration via Ehrenfest’s theorem [84], i.e.,

[TABLE]

In Fig. 3 we show the computed high-order harmonic spectrum computed for a hydrogen atom under a laser pulse with carrier frequency corresponding to $800\text{\,}\mathrm{nm}$, intensity full width at half maximum $4.8\text{\,}\mathrm{fs}$, and electric field peak intensity $36\text{\,}\mathrm{GV}\text{/}\mathrm{m}$. The spectrum is a typical HHG spectrum from a femtosecond laser pulse, and it ends at the cut-off energy $\sim$ $2\text{\,}\mathrm{a.u.}$ in a decent agreement with the three-step model’s cutoff $\sim$ $1.7\text{\,}\mathrm{a.u.}$ [80, 81].

4.3 Metal nanotips and inhomogeneous fields

Metal nanotips have recently attracted attention due to their ability to enhance the laser electric field at the tip apex via plasmonic effects (see, e.g., [19, 20, 21, 22, 23, 24, 25, 26, 27, 28]). Fiend is capable of computing single-electron dynamics of these metal nanotapers including correctly the interaction with the inhomogeneous plasmon enhanced laser field.

The first step is to compute the plasmon enhanced near field. Here we consider a gold nanotip with apex radius $7\text{\,}\mathrm{nm}$ and full opening angle of $20\,{}^{\circ}$ as demonstrated in Fig. 4(a). For typical laser wavelengths, such as $800\text{\,}\mathrm{nm}$ used here, a quasistatic description of the laser vector potential is applicable, i.e.,

[TABLE]

where the spatial form $\mathbf{A}_{s}$ is the same as for the electric field

[TABLE]

The electrostatic potential $U(\mathbf{r})$ can be computed from the Poisson equation

[TABLE]

where $\epsilon(\mathbf{r})$ is the dielectric function of the material at position $\mathbf{r}$. Inside the gold nanotip the complex dielectric function (at $800\text{\,}\mathrm{n}\mathrm{m}$) is $\epsilon_{\mathrm{Au}}=-24.061+1.5068\,{\mathrm{i}\mkern 1.0mu}$ [85, 86] and at vacuum $\epsilon_{\mathrm{vac}}=1$.

We can describe Eq. (31) using cylindrical coordinates. The interface condition

[TABLE]

gives rise to the weak form

[TABLE]

with the test function space

[TABLE]

and the trial function space

[TABLE]

The simulation domain is a rectangular area, and the ZDBCs described above are imposed on the external boundaries denoted by $\Gamma_{D}$. Note that Fig. 4(a) demonstrates the domain and boundaries for the quantum simulation – for Poisson equation we employ a much larger simulation domain. The equation (31) can be solved easily with the FEniCS FEM-suite and we provide an example script at demos/nanotip/1_near_field.py

The computed spatial distribution of the field, $\|A_{s}(x,y,z)\|$, is demonstrated in Fig. 4(b). The field is at its maximum at the nanotip apex and decays rapidly with the distance to the tip. We note that this quasistatic solution provides field-enhancement factor of $f\approx 10$ which is comparable with experimental results [87].

In TDSE, we use a potential well for the static potential,

[TABLE]

where $\Phi=$5.31\text{,}\mathrm{e}\mathrm{V}$$ is the work function of Au (111)-surface [88].

We prepare the electron as a Gaussian wave packet [Fig. 4(c)] and assemble the system matrices – including the inhomogeneous field – with demos/nanotip/3_prepare_tdse.py. The preparation step is slightly more complex than for linearly polarized vector potentials, but the provided example script should work as a template for further expansions.

We use a $800\text{\,}\mathrm{nm}$ pulse with $15\text{\,}\mathrm{fs}$ full width at half maximum for the Gaussian envelope. The maximum electric field amplitude without nanostructure enhancement is set to $30\text{\,}\mathrm{G}\mathrm{V}\mathrm{/}\mathrm{n}\mathrm{m}$. Snapshots of the electron density during propagation are shown in Figs. 4(c)-(f). As expected, the electron emission is concentrated at the tip apex where the field enhancement is at its highest. Only a small percentage (\sim$$2 %) of the electron density is absorbed at the simulation box boundary.

5 Summary

We have presented Fiend – a versatile solver for single-particle quantum dynamics in cylindrically symmetric systems. This Python package provides an easy path for the study of nonlinear strong-field phenomena in atoms and nanostructures under both homogeneous and inhomogeneous vector potentials. Moreover, the finite element discretization used by Fiend can adapt to complicated system geometries. The package is parallelized using MPI, and much of the high-performance computing is delegated to state-of-the-art numerical libraries. Fiend is modular and easy to extend, and we provide comprehensive documentation for the program code.

We have demonstrated the capabilities of Fiend by simulating two different scenarios. The simulation of high-order harmonic generation from a laser-driven hydrogen atom shows that Fiend can be applied to traditional strong-field phenomena in atoms. Furthermore, a simulation of the photoionization of a gold nanotip demonstrates Fiend’s suitability for studying plasmon-enriched strong-field phenomena in nanostructures.

Acknowledgements

We are grateful to Matti Molkkari and Joonas Keski-Rahkonen for insightful discussions. This work was supported by the Alfred Kordelin foundation and the Academy of Finland (grant no. 304458). We also acknowledge CSC – the Finnish IT Center for Science – for computational resources.