Finite-difference time-domain simulation of strong-field ionization: Perfectly matched layer approach

TL;DR

This paper introduces an FDTD simulation method with PMLs for modeling strong-field electron ionization, demonstrating significant resource savings and superior performance over existing techniques for short-range potentials.

Contribution

It develops a PML-based FDTD approach for the time-dependent Schrödinger equation, showing improved efficiency and accuracy in simulating electron ionization under strong fields.

Findings

PMLs reduce computational domain size by several orders of magnitude.

PMLs outperform ECS for short-range potentials in FDTD simulations.

The method's accuracy is validated against known numerical and analytical results.

Abstract

A Finite-Difference Time-Domain (FDTD) scheme with Perfectly Matched Layers (PMLs) is considered for solving the time-dependent Schr\"{o}dinger equation, and simulate the ionization of an electron initially bound to a one-dimensional $\delta$-potential, when applying a strong time-oscillating electric field. The performance of PMLs based on different absorption functions are compared, where we find slowly growing functions to be preferable. PMLs are shown to be able to reduce the computational domain, and thus the required numerical resources, by several orders of magnitude. This is demonstrated by testing the proposed method against an FDTD approach without PMLs and a very large computational domain. We further show that PMLs outperform the well known Exterior Complex Scaling (ECS) technique for short-range potentials when implemented in FDTD, though ECS remains superior for long-range potentials. The accuracy of the method is furthermore demonstrated by comparing with known numerical and analytical results for the $\delta$-potential.

Equations

$$\displaystyle i\frac{d}{dt}\psi\mathopen{}\mathclose{{}\left(\boldsymbol{r},t}\right)=\mathopen{}\mathclose{{}\left[-\frac{1}{2}\nabla^{2}+H_{g}\mathopen{}\mathclose{{}\left(\boldsymbol{r},t}\right)+V\mathopen{}\mathclose{{}\left(\boldsymbol{r}}\right)}\right]\psi\mathopen{}\mathclose{{}\left(\boldsymbol{r},t}\right)\thinspace,$$

$$\displaystyle H_{V}=\boldsymbol{p}\cdot\boldsymbol{A}\mathopen{}\mathclose{{}\left(t}\right)\thinspace,$$

$$\displaystyle\boldsymbol{A}\mathopen{}\mathclose{{}\left(t}\right)=A_{0}\cos\mathopen{}\mathclose{{}\left(\omega t}\right)\hat{\boldsymbol{x}}\thinspace.$$

$$\displaystyle\mathbcal{E}\mathopen{}\mathclose{{}\left(t}\right)=-\frac{\partial\boldsymbol{A}\mathopen{}\mathclose{{}\left(t}\right)}{\partial t}=\mathcal{E}_{0}\sin\mathopen{}\mathclose{{}\left(\omega t}\right)\hat{\boldsymbol{x}}$$

$$\displaystyle H_{L}=\boldsymbol{r}\cdot\mathbcal{E}\mathopen{}\mathclose{{}\left(t}\right)\thinspace.$$

$$\displaystyle\psi\mathopen{}\mathclose{{}\left(\boldsymbol{r},t}\right)=\psi_{\mathrm{ex}}\mathopen{}\mathclose{{}\left(\boldsymbol{r},t}\right)\quad\text{for }\mathopen{}\mathclose{{}\left|\boldsymbol{r}}\right|\leq R_{0}\thinspace,$$

$$\displaystyle\sigma\mathopen{}\mathclose{{}\left(R_{0}}\right)=1-\frac{\mathopen{}\mathclose{{}\left|\mathopen{}\mathclose{{}\left<\psi_{\mathrm{ex}}|\psi}\right>_{R_{0}}}\right|^{2}}{\mathopen{}\mathclose{{}\left<\psi_{\mathrm{ex}}|\psi_{\mathrm{ex}}}\right>_{R_{0}}\mathopen{}\mathclose{{}\left<\psi|\psi}\right>_{R_{0}}}\thinspace,$$

$$\displaystyle\mathopen{}\mathclose{{}\left<f|g}\right>_{R_{0}}=\int_{\mathopen{}\mathclose{{}\left|\boldsymbol{r}}\right|\leq R_{0}}f^{*}\mathopen{}\mathclose{{}\left(\boldsymbol{r}}\right)g\mathopen{}\mathclose{{}\left(\boldsymbol{r}}\right)d\boldsymbol{r}\thinspace.$$

$$\displaystyle x\rightarrow\tilde{x}=\begin{cases}x\quad&\mathrm{for}\,\mathopen{}\mathclose{{}\left|x}\right|<R_{0}\\ e^{i\theta}\mathopen{}\mathclose{{}\left(x\pm R_{0}}\right)\mp R_{0}\quad&\mathrm{for}\,\mp x>R_{0}\thinspace,\end{cases}$$

$$\exp\mathopen{}\mathclose{{}\left(-iH_{L}\Delta t}\right)=\exp\mathopen{}\mathclose{{}\left\{-i\mathcal{E}\mathopen{}\mathclose{{}\left(t}\right)\mathopen{}\mathclose{{}\left[\cos\theta\mathopen{}\mathclose{{}\left(x\pm R_{0}}\right)\mp R_{0}}\right]\Delta t}\right\}\\ \times\exp\mathopen{}\mathclose{{}\left\{\mathcal{E}\mathopen{}\mathclose{{}\left(t}\right)\sin\theta\mathopen{}\mathclose{{}\left(x\pm R_{0}}\right)\Delta t}\right\}\thinspace,$$

$$\exp\mathopen{}\mathclose{{}\left(-iH_{V}\Delta t}\right)=\exp\mathopen{}\mathclose{{}\left(-i\cos\theta A\mathopen{}\mathclose{{}\left(t}\right)p_{x}\Delta t}\right)\\ \times\exp\mathopen{}\mathclose{{}\left(-\sin\theta A\mathopen{}\mathclose{{}\left(t}\right)p_{x}\Delta t}\right)\thinspace,$$

$$\displaystyle x\to\tilde{x}=\begin{cases}x\quad&\mathrm{for}\,\mathopen{}\mathclose{{}\left|x}\right|<R_{0}\\ x+i\sigma_{0}\int^{x}f\mathopen{}\mathclose{{}\left(x^{\prime}}\right)dx^{\prime}\quad&\mathrm{for}\,\mathopen{}\mathclose{{}\left[x}\right]>R_{0}\thinspace,\end{cases}$$

$$\displaystyle\frac{\partial}{\partial x}\to c\mathopen{}\mathclose{{}\left(x}\right)\frac{\partial}{\partial x}\thinspace,$$

$$\displaystyle i\frac{\partial}{\partial t}\psi=\mathopen{}\mathclose{{}\left[-\frac{1}{2}c\mathopen{}\mathclose{{}\left(x}\right)\frac{\partial}{\partial x}c\mathopen{}\mathclose{{}\left(x}\right)\frac{\partial}{\partial x}+H_{g}\mathopen{}\mathclose{{}\left(x,t}\right)+V\mathopen{}\mathclose{{}\left(x}\right)}\right]\psi\thinspace,$$

$$\displaystyle f\mathopen{}\mathclose{{}\left(y}\right)=\Theta\mathopen{}\mathclose{{}\left(y}\right)\begin{cases}d/\mathopen{}\mathclose{{}\left(d-y+\epsilon}\right)-1\\ y^{2}\\ y^{3}\\ \tanh\mathopen{}\mathclose{{}\left(2y/d-1}\right)-\tanh\mathopen{}\mathclose{{}\left(-1}\right)\thinspace,\end{cases}$$

$$\displaystyle V\mathopen{}\mathclose{{}\left(x}\right)=\begin{cases}-\frac{1}{2b}\quad&\mathrm{for}\,\mathopen{}\mathclose{{}\left|x}\right|<b\thinspace,\\ 0\quad&\mathrm{otherwise}\thinspace.\end{cases}$$

$$\displaystyle\alpha\mathopen{}\mathclose{{}\left(\omega}\right)=\frac{2-\omega^{2}-\sqrt{1+2\omega}-\sqrt{1-2\omega}}{\omega^{4}}\thinspace.$$

$$\displaystyle P_{\mathrm{bound}}\mathopen{}\mathclose{{}\left(t;\mathcal{E}_{0}}\right)=\sum_{b}\mathopen{}\mathclose{{}\left|\mathopen{}\mathclose{{}\left<\varphi_{b}|\psi_{\mathrm{ex}}}\right>}\right|^{2}\thinspace,$$

$$\displaystyle\mathopen{}\mathclose{{}\left|\varphi_{B}\mathopen{}\mathclose{{}\left(x}\right)}\right|\geq\mathopen{}\mathclose{{}\left|\varphi_{b}\mathopen{}\mathclose{{}\left(x}\right)}\right|\quad\text{for }\mathopen{}\mathclose{{}\left|x}\right|\geq L\thinspace,$$

$$\displaystyle\Gamma\mathopen{}\mathclose{{}\left(t;\mathcal{E}_{0}}\right)=-\frac{d}{dt}\ln P_{\mathrm{bound}}\mathopen{}\mathclose{{}\left(t;\mathcal{E}_{0}}\right)\thinspace.$$

$$\displaystyle\mathopen{}\mathclose{{}\left<\Gamma\mathopen{}\mathclose{{}\left(\mathcal{E}_{0}}\right)}\right>=\frac{\omega}{2n\pi}\int_{t_{0}}^{t_{0}+2n\pi/\omega}\Gamma\mathopen{}\mathclose{{}\left(t;\mathcal{E}_{0}}\right)dt\thinspace,$$

$$\displaystyle\frac{\mathcal{E}_{\mathrm{DC}}^{1/3}}{2^{2/3}\pi}-\text{Ai}\mathopen{}\mathclose{{}\left(\lambda}\right)\text{Bi}\mathopen{}\mathclose{{}\left(\lambda}\right)=i\text{Ai}^{2}\mathopen{}\mathclose{{}\left(\lambda}\right)\thinspace,$$

$$\displaystyle\mathopen{}\mathclose{{}\left<\Gamma_{\mathrm{Adiabatic}}\mathopen{}\mathclose{{}\left(\mathcal{E}_{0}}\right)}\right>=\frac{\omega}{2\pi}\int_{0}^{2\pi/\omega}\Gamma_{\mathrm{DC}}\mathopen{}\mathclose{{}\left[\mathcal{E}\mathopen{}\mathclose{{}\left(t}\right)}\right]dt\thinspace.$$

$$\displaystyle\Gamma_{\mathrm{asymp.}}\mathopen{}\mathclose{{}\left(\mathcal{E}_{\mathrm{DC}}}\right)=\mathopen{}\mathclose{{}\left(1-\frac{5}{3}\mathcal{E}_{\mathrm{DC}}}\right)\exp\mathopen{}\mathclose{{}\left(-\frac{2}{3\mathcal{E}_{\mathrm{DC}}}}\right)\thinspace.$$

$$\displaystyle\langle\Gamma_{\mathrm{PPT}}\mathopen{}\mathclose{{}\left(\mathcal{E}_{0}}\right)\rangle=\mathopen{}\mathclose{{}\left(\frac{3\mathcal{E}_{0}}{\pi}}\right)^{1/2}\exp\mathopen{}\mathclose{{}\left(-\frac{2}{3\mathcal{E}_{0}}}\right)\thinspace.$$

$$\displaystyle V\mathopen{}\mathclose{{}\left(x}\right)=-\frac{1}{\sqrt{x^{2}+2}}\thinspace.$$

$$i\frac{\psi\mathopen{}\mathclose{{}\left(x,t_{j+1}}\right)-\psi\mathopen{}\mathclose{{}\left(x,t_{j}}\right)}{\Delta t}\\ =\frac{1}{2}\mathopen{}\mathclose{{}\left[H\mathopen{}\mathclose{{}\left(t_{j+1}}\right)\psi\mathopen{}\mathclose{{}\left(x,t_{j+1}}\right)+H\mathopen{}\mathclose{{}\left(t_{j}}\right)\psi\mathopen{}\mathclose{{}\left(x,t_{j}}\right)}\right]\thinspace,$$

$$T_{\mathrm{PML}}=-\frac{1}{2}c\mathopen{}\mathclose{{}\left(x}\right)\frac{\partial}{\partial x}c\mathopen{}\mathclose{{}\left(x}\right)\frac{\partial}{\partial x}\psi\\ =-\frac{1}{2}\mathopen{}\mathclose{{}\left[c^{2}\mathopen{}\mathclose{{}\left(x}\right)\frac{\partial^{2}\psi}{\partial x^{2}}+c\mathopen{}\mathclose{{}\left(x}\right)\frac{\partial c\mathopen{}\mathclose{{}\left(x}\right)}{\partial x}\frac{\partial\psi}{\partial x}}\right]\thinspace.$$

$$\frac{\partial \psi }{\partial x}$$

$$\frac{{\partial }^2\psi }{\partial x^2}$$

$$\displaystyle\psi\mathopen{}\mathclose{{}\left(\tilde{x}}\right)\approx\sum_{n=-p}^{p}l_{n}\mathopen{}\mathclose{{}\left(\tilde{x}}\right)\psi\mathopen{}\mathclose{{}\left(\tilde{x}_{n}}\right)$$

Full text

I Introduction

It has long been established that excitons must be taken into account to accurately describe the optical properties of solids. If the excitons are strongly bound, further complications arise as the excitons must be dissociated into free charge carriers before they can supply an electrical current. Interest in using monolayer transition metal dichalcogenides (TMDs) in electronic devices has increased dramatically during recent years. These materials are known to support strongly bound excitons Wang et al. (2012); Geim and Grigorieva (2013); Olsen et al. (2016), which dominate their linear and nonlinear optical properties Wang et al. (2012); Ramasubramaniam (2012); Qiu et al. (2013); Trolle et al. (2015). However, efficient generation of photocurrents in e.g. photodetectors and solar cells require dissociation of excitons into free electrons and holes. It is therefore of significant interest to obtain efficient methods of inducing dissociation of bound excitons in materials such as TMDs. Recently, static in-plane electric fields have proven a promising candidate to aid dissociation Massicotte et al. (2018); Haastrup et al. (2016); Pedersen et al. (2016a); Kamban and Pedersen , and the rates induce by such fields are readily calculated by the well known complex scaling method Balslev and Combes (1971); Aguilar and Combes (1971); Haastrup et al. (2016); Kamban and Pedersen , or by using a hypergeometric resummation technique Mera et al. (2015). The situation is not as simple for time-dependent fields, however, which motivates the search for efficient methods of calculating the dissociation rates induced by alternating fields. The problem is remarkably similar to calculating the ionization rate of atoms Perelomov et al. (1966); Keldysh (1965); Gersten and Mittleman (1974); Tolstikhin et al. (2011); Trinh et al. (2013); Greenwood and Eberly (1991); Javanainen et al. (1988); Plummer et al. (1998); Tang and Shakeshaft (1987); Dörr and Shakeshaft (1988), for which Floquet theory Floquet (1883) implemented with complex scaling Holt et al. (1983); Maquet et al. (1983) has proven very useful. A major drawback with traditional Floquet theory, however, is that it only applies to periodic fields. In the present paper, we seek to develop a different approach based on PMLs. To accurately described TMD excitons, one must use the Keldysh potential Keldysh (1979); Trolle et al. (2017). As this potential is rather complicated, we shall test the method developed here by calculating the ionization rate of an electron bound to the one-dimensional zero-range potential -\delta\mathopen{}\mathclose{{}\left(x}\right).

One of the most reliable techniques of obtaining accurate results in intense laser-matter interactions is to propagate the time-dependent Schrödinger equation (TDSE) and calculate the relevant observables. Strong electric fields lead to a large probability flux traveling out of the central region occupied by the initially localized wave function, which, of course, is exactly what we mean by ionization. If one simply uses Dirichlet boundary conditions at the edges of the simulation domain, a huge domain is required to avoid reflections of the wave function from the boundary. Several methods have been designed to circumvent this problem, with one of the most common ones being absorbing boundaries outside a specified interior domain. Popular methods include complex absorbing potentials (CAPs) Riss and Meyer (1996), absorbing masks Krause et al. (1992), and ECS He et al. (2007); Tao et al. (2009). The goal of the absorbing layers is to leave the wave function unaltered in the interior region, while absorbing it as it moves out of this region. If it is absorbed sufficiently quickly, then one is able to use Dirichlet boundary conditions at a distance into the layer with minimal flux reaching this point and thus avoiding spurious reflections.

The ECS method has been given a lot of attention in recent years, and rightfully so, as it has been shown to be a very efficient absorber in time-dependent Schrödinger problems McCurdy et al. (1991); Scrinzi (2010). A related method that has been given less attention in quantum mechanics is the use of PMLs. The PML method was developed by Berenger for solving Maxwell’s equations Berenger (1994), and has since been used extensively in classical electromagnetism, where PMLs are applied efficiently in FDTD Inan and Marshall (2011); Hagness and Taflove (2005), in frequency-domain finite-element Jin (2002), and in Fourier-series Zhang et al. (2008) approaches. Lu and Zhu additionally proposed a perturbative approach Lu and Zhu (2005) to deal with undesired effects of the PML when simulating optical wave guides, and PMLs have, furthermore, been utilized to study sound waves Zuo and Fan (2017). Given the success of the PML method in solving problems in electromagnetism, interest in applying it to Schrödinger problems has slowly been increasing. Zheng used it to solve the nonlinear Schroödinger equation Zheng (2007), and Nissen and Kreiss have since tried to optimize the PML method for the Schrödinger equation with time-independent potentials Nissen and Kreiss (2011), and have, together with Karlsson, applied it to a reactive scattering problem Nissen et al. (2010). PMLs have also been applied to time-dependent-density-functional theory (TDDFT) Lehtovaara et al. (2011), and the Dirac equation Pinaud (2015). It is, however, surprising how comparatively little work has been done on applying PMLs in Schrödinger problems, in particular, for explicitly time-dependent problems, such as intense laser-matter interactions.

In the present paper, we develop a method based on a finite-difference time-domain scheme including a PML (FDTD-PML) to describe the ionization of an electron bound by the zero-range $\delta$-potential. This potential has previously been used to, e.g., study the optical response in one-dimensional semiconductors Pedersen (2015) and to model ionization of the $H^{-}$ ion Scharf et al. (1991). The convergence of the method will be compared to a standard FDTD scheme using Dirichlet boundary conditions, as well the well known ECS method. We will show that for a short-range potential, the PML method far outperforms ECS when both are implemented as finite-difference schemes. Subsequently, the method will be used to analyze limiting cases, where analytical results can obtained Fernández and Castro (1985); Maize and Williams (2004); Postma (1984). This is done in order to check that the FDTD-PML results remain physical, even though the wave function will be absorbed by the PML. The ionization rate is thereafter calculated as a function of frequency and field strength. Finally, we show that PMLs are not well suited for potentials that reach far into the absorbing layers.

II Electron in a Laser Field

We seek to solve the time-dependent Schrödinger equation for an electron, initially bound to a localized potential, perturbed by a monochromatic laser field (atomic units are used throughout)

[TABLE]

where $H_{g}$ describes the interaction between the electron and the laser field. Here, the subscript $g$ refers to the gauge, in which the interaction is considered. We shall work only in the dipole approximation such that neither electric fields nor vector potentials have any spatial dependence. This leads to the interaction in the velocity gauge (VG) being

[TABLE]

where $\boldsymbol{p}$ is the momentum operator and $\boldsymbol{A}$ is the vector potential. Note that the usual A\mathopen{}\mathclose{{}\left(t}\right)^{2}/2 term has been removed by a unitary transformation. We will consider the monochromatic field defined by

[TABLE]

In the length gauge (LG), the interaction is given in terms of the electric field

[TABLE]

by

[TABLE]

The goal is to be able to reproduce the exact wave function in an interior box \mathopen{}\mathclose{{}\left|\boldsymbol{r}}\right|\leq R_{0} for relevant time periods. That is, we seek to modify the TDSE so that the solution to the modified equation $\psi$ satisfies

[TABLE]

where $\psi_{\mathrm{ex}}$ is the exact wave function. To be able to quantify the error by a single number, we will use the error measurement introduced by Scrinzi Scrinzi (2010)

[TABLE]

where the scalar product is to be taken in the region \mathopen{}\mathclose{{}\left|\boldsymbol{r}}\right|\leq R_{0}, i.e.

[TABLE]

III Exterior Complex Scaling

The literature covering complex scaling is vast (see e.g. Balslev and Combes (1971); Aguilar and Combes (1971); Doolen et al. (1974); Ho (1981); Reed and Simon (1982); Bengtsson et al. (2008) and references therein), and for this reason we shall only describe briefly the most relevant aspects to the present paper before describing how PMLs are implemented. For simplicity, we will restrict the discussion to one dimension $x$. The simplest form of complex scaling is implemented by scaling the coordinates uniformly according to $x\rightarrow xe^{i\theta}$ , where $\theta$ will be taken as a purely real number. The motivation is that the outgoing waves \exp\mathopen{}\mathclose{{}\left(ikx}\right) become exponentially decaying waves if the rotational angle $\theta$ is chosen large enough. This transformation, referred to as uniform complex scaling (UCS), has been used with great success in finding ionization rates for static electric fields Herbst and Simon (1978); Pedersen et al. (2016b); Massicotte et al. (2018) and in solving the TDSE Bengtsson et al. (2008). In certain situations, however, one may wish to leave the domain untransformed in an interior region and introduce complex scaling only in an outer region. The original motivation was that UCS cannot be used with the Born-Oppenheimer approximation Simon (1979). Furthermore, when dealing with sufficiently weak electric fields, ionization rates typically become so low that UCS results in a wave function that (numerically) vanishes before it reaches the important region far from the core Kamban and Pedersen . The ECS Simon (1979); McCurdy et al. (2004, 1991); Rescigno and McCurdy (2000) procedure circumvents these problems, and is implemented in one dimension by the transformation

[TABLE]

where $R_{0}$ is referred to as the scaling radius. This transformation turns outgoing waves into decaying waves in the absorbing layer while leaving them unaffected in the interior region.

The resulting behaviour of the wave function in the absorbing layer is slightly different in the two gauges. The exponential propagator can be constructed as usual (see Ref. He et al. (2007)), and the perturbing part can be written as

[TABLE]

in LG and as

[TABLE]

in VG. The first terms on the right-hand sides of Eqs. (10) and (11) are oscillatory, while the second terms are either exponentially increasing or exponentially decreasing. In LG, this depends on the sign of the oscillatory field \mathcal{E}\mathopen{}\mathclose{{}\left(t}\right) and in VG on the sign of A\mathopen{}\mathclose{{}\left(t}\right). Thus, in both cases one may obtain an undesired exponentially increasing behavior in the absorbing layer. Given that the behavior depends on the sign of the field or vector potential, the propagators will oscillate between amplifying and damping the wave function exponentially. In practice, we have found that the exponential behavior outside of the scaling radius is much more apparent in LG than inVG. Furthermore, in LG, the exponential behavior is more volatile for larger $x$, which may lead to numerical instabilities if a wide absorbing layer is desired. In practice, we have not found these growing terms to cause numerical instabilities for moderate frequencies, while for low frequencies they lead to a numerically diverging wave function. This is in agreement with the observations in Ref. He et al. (2007).

IV Perfectly Matched Layers

The PML scheme for the TDSE is usually derived by assuming that the potential is both spatially and temporally invariant, and then modal analysis is performed on the Laplace-transformed equation to ensure that the solution decays outside the interior domain, i.e. \mathopen{}\mathclose{{}\left|x}\right|>R_{0} Zheng (2007); Nissen and Kreiss (2011); Nissen et al. (2010). The transformation can be formulated as

[TABLE]

where $\sigma_{0}$ is a constant referred to as the absorption strength and $f$ is the absorption function. The absorption function is zero inside the interior \mathopen{}\mathclose{{}\left|x}\right|\leq R_{0} and positive otherwise. Specific forms will be discussed later. Unlike in ECS, the transformation in Eq. 12 is not applied to the potential. Thus, the PML method can be understood as a transformation of the differential operator

[TABLE]

where c\mathopen{}\mathclose{{}\left(x}\right)=1/\mathopen{}\mathclose{{}\left[1+i\sigma_{0}f\mathopen{}\mathclose{{}\left(x}\right)}\right] . The PML equation in one dimension therefore becomes

[TABLE]

which coincides with the usual TDSE inside a box of radius $R_{0}$ as c\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left|x}\right|\leq R_{0}}\right)=1. Note that the momentum operator in $H_{V}$ (see Eq. 2) is also transformed according to Eq. 13. As the transformation is only applied to the spatial derivatives, it is only reasonable to expect Eq. (14) to yield a good approximation if $R_{0}$ is chosen sufficiently large so that the variations in the potential V\mathopen{}\mathclose{{}\left(x}\right) in the exterior are negligible. This is the case for any non-zero $R_{0}$ for the zero-range potential V\mathopen{}\mathclose{{}\left(x}\right)=-\delta\mathopen{}\mathclose{{}\left(x}\right). However, the interaction in LG effectively modifies the potential so that it includes a linear term which does not vanish outside the interior. One may therefore speculate to which degree Eq. 14 in LG is able to approximate the exact wave function in the interior. Indeed, as we show numerically below, implementing the PML in LG introduces significantly larger errors than in VG. The technical details of implementing both the ECS and PML method can be found in App. A.

IV.1 Absorption function

It is important that the absorption function $f$ be chosen positive to ensure decay of the wave function as it travels out of the interior domain. Previous choices include low-degree power functions Zheng (2007); Berenger (1994) and singular functions Bermúdez et al. (2004). It is interesting to examine whether or not there is a substantial difference between the numerical accuracy obtainable with the different functions. To this end, we will compare four different absorption functions, namely

[TABLE]

where y=\mathopen{}\mathclose{{}\left|x}\right|-R_{0}, $\epsilon$ is some small positive number, and $\Theta$ is a step function equal to unity for $y\geq 0$ and zero otherwise, and $d$ is the width of the absorbing layer. For $\epsilon$, we have used $10^{-4}$, as we have not found the results to be highly dependent on $\epsilon$.

V Short-range Potential

The first potential we will consider is the short-range potential

[TABLE]

This potential can be seen as a discrete approximation to the zero-range potential -\delta\mathopen{}\mathclose{{}\left(x}\right), for which we want to calculate the ionization rate. We have used $b=5\times 10^{-3}$ for the calculations in the present paper, which leads to a ground-state energy of $E_{0}=-0.4967\,\mathrm{a.u.}$ (as opposed to $-1/2\,\mathrm{a.u.}$ for the zero-range potential).

Figure 1 shows the absolute square of the wave function calculated with the ECS and the PML method at three different times. The oscillations in LG outside the scaling radius ($R_{0}=20\,\mathrm{a.u.}$) that were discussed briefly in the ECS section above are immediately clear. The absolute square of the wave function is, however, graphically indistinguishable in the interior for ECS implemented in LG and VG. For the PML calculations, the results are not as equal-footed. The LG calculation introduces non-negligible reflections leading to a large error inside the interior. This is, of course, what we seek to avoid and thus one must be careful in implementing PML in LG. For VG, the PML and ECS wave functions are indistinguishable for \mathopen{}\mathclose{{}\left|x}\right|\leq R_{0}.

To perform error analysis, we need a reference function. It was obtained by calculating the error in Eq. 7 inside $R_{0}=20$ a.u. at $t=200$ a.u. between two wave functions without any transformation and with Dirichlet boundary conditions set at $x=\pm n\times 500$ for $\psi$ and x=\pm\mathopen{}\mathclose{{}\left(n+1}\right)\times 500 a.u. for $\psi_{\mathrm{ex}}$, where $n$ is a positive integer, which was increased by one until the error vanished within numerical precision. This occurred at $n=8$, and thus without any absorbing layer a domain width of at least $8000$ a.u. is needed. To ensure a numerically exact reference function, we have used a domain width of $10000$ a.u. for $\psi_{\mathrm{ex}}$ in the error calculations. In Fig. 2 we show the error calculated by Eq. 7 at $t=200$ a.u. as a function of the PML width $d$ for the different absorption functions with various absorption strengths $\sigma_{0}$. It is clear that the PML should be implemented in VG to obtain an accurate wave function in the interior. We also notice that an absorption function that grows slowly leads to a lower error at the cost of slower convergence. The reason is that less of the wave function will be (numerically) reflected upon entering the absorbing layer when the transition is more gradual. The error introduced by the power functions converge to a constant value for a sufficiently large $d$. This indicates that the entire outgoing flux has either been absorbed or reflected from the layer, and further increasing $d$ does not make any difference. It is worth noting that the PML method with a square absorption function (Fig. 2(b)) leads to errors of order $10^{-15}$ for an absorbing layer of around $d=40$ a.u.. Thus, the domain width required to obtain an excellent approximation to the exact wave function is 2\mathopen{}\mathclose{{}\left(d+R_{0}}\right)=120 a.u., significantly lower than the domain width of $8000$ a.u. needed without any absorbing layer. For the nearly singular function and the $\tanh$ function, increasing $d$ leads to an absorption function that grows more slowly. For this reason, the errors they induce do not converge in the same manner as those induced by the power functions. Rather, they continue to decrease as the absorbing layers become wider.

The PML errors should be compared to those introduced by the ECS method, shown in Fig. 3, for which both LG and VG calculations converge to the same error, and are nearly indistinguishable. A low rotational angle $\theta$ can be seen to introduce lower errors, as less of the wave function will be reflected upon entering the absorbing layer, exactly as with the PML method. By comparing the errors introduced by ECS to those for PMLs with a square absorption function in Fig. 2(b), we see that the ECS errors for comparable absorption widths converge to values that are around six orders of magnitude larger. For PMLs with a quadratic absorption function (Fig. 2(c)), the ECS errors are around four orders of magnitude larger. This might be due to the poor performance of ECS when implemented in finite-difference schemes McCurdy et al. (2004). The same behavior is observed in Fig. 4, where the error is shown as a function of time. The PML LG calculation introduces much larger errors than the other two methods, and the ECS LG and VG errors are graphically indistinguishable. Again, the PML VG calculation leads to the lowest error by several orders of magnitude. The dotted lines show a scaling radius of $R_{0}=10\,\mathrm{a.u.}$ as opposed to $R_{0}=20\,\mathrm{a.u.}$. As can be seen, reducing the size of the box does not have a large impact on the errors in the interior domain.

V.1 Polarizability and ionization

In the previous section, the error of the wave function inside the box \mathopen{}\mathclose{{}\left|x}\right|<R_{0} was analyzed. If the wave function can be reproduced, then the desirable observables can be calculated, as long as the box size is chosen adequately. While the error measurement defined by Eq. 7 is a meaningful parameter, we are unable to directly relate it to physical observables. As an additional check, we therefore demonstrate that we are able to reproduce the frequency dependent polarizability in the weak-field limit. Numerically, the polarizability $\alpha$ can be found by calculating \mathopen{}\mathclose{{}\left<x}\right>/\mathcal{E}_{0} in the weak field limit, and relating it to the real and imaginary part of $\alpha$. Here, \mathopen{}\mathclose{{}\left<x}\right> is the average value of $x$. As the field amplitude is extremely low, and \mathopen{}\mathclose{{}\left<x}\right> is only needed over a single period, the integral in \mathopen{}\mathclose{{}\left<x}\right> can, to an excellent approximation, be restricted to the interior region. An analytical expression can be found for $\alpha$ of the $\delta$-potential ground-state using linear perturbation theory Postma (1984). It is given by

[TABLE]

As can be seen in Fig. 5, the PML simulations using a low field strength of $\mathcal{E}_{0}=10^{-6}$ a.u. are in excellent agreement with the analytical results.

A special interest in the present paper is to obtain the strong-field ionization rate. The probability of occupying a bound state (that is, not being ionized) can be found by

[TABLE]

where the sum is to be taken over all bound states. For numerical calculations, however, one may cut the sum after convergence to a desired number of digits. Let us denote the most delocalized state included in the sum by $\varphi_{B}$, such that from some number $L$

[TABLE]

where $b$ refers to all states included. If $L$ is chosen such that $\varphi_{B}$ is negligible for \mathopen{}\mathclose{{}\left|x}\right|>L, then all integrals in Eq. 18 may be restricted to \mathopen{}\mathclose{{}\left|x}\right|<L. By choosing the scaling radius $R_{0}$ to coincide with $L$, we can therefore describe all bound states, as well as obtain an excellent approximation to the wave function, in the interior domain, allowing us to implement Eq. 18 in the present approach. For the short-range potential defined by Eq. 16 there is only one bound state and it decays exponentially for \mathopen{}\mathclose{{}\left|x}\right|>b. Therefore, a scaling radius of $R_{0}=20\,\mathrm{a.u.}$ is expected to be adequate. The probability of not being ionized can be seen as a function of time in Fig. 6, where the ECS and PML calculations are compared to a converged Crank-Nicolson calculation in an untransformed domain. As is evident, the result obtained by the PML method in LG is the only one that can be distinguished from the other ones. This is yet another indication that care must be taken when implementing PMLs in LG.

To obtain the time-dependent ionization rate $\Gamma$, we use

[TABLE]

The ionization rate defined by Eq. 20 will oscillate in time. It is therefore convenient to average the time-dependent ionization rate over a number of periods to remove these oscillations, and thereby obtain a time-independent ionization rate, i.e.

[TABLE]

where $n$ is the number of periods, and $t_{0}$ is an initial time taken after the field has been turned on. The ionization rate averaged over two periods can be seen in the upper panel of Fig. 7 for $\omega=0.2\,\mathrm{a.u.}$. As is evident, the four methods yield identical results. Furthermore, the shape of the ionization rate is consistent with the results for three-photon ionization in Ref. Scharf et al. (1991).

In the adiabatic limit, one can obtain analytical results for the ionization rate of the zero-range potential. By setting up the Schrödinger equation for a static electric field $\mathcal{E}_{\mathrm{DC}}>0$ and requiring that the wave function becomes an outgoing wave as $x\rightarrow-\infty$ Fernández and Castro (1985), one can obtain the following condition

[TABLE]

where $\lambda=-2^{1/3}E\mathcal{E}_{\mathrm{DC}}^{-2/3}$, and Ai and Bi are Airy functions of the first and second kind Abramowitz and Stegun (1972), respectively. Solving Eq. 22 numerically one obtains complex energies and the DC ionization rate is then given by \Gamma_{\mathrm{DC}}\mathopen{}\mathclose{{}\left(\mathcal{E}_{\mathrm{DC}}}\right)=-2\text{Im}\mathopen{}\mathclose{{}\left[E\mathopen{}\mathclose{{}\left(\mathcal{E}_{\mathrm{DC}}}\right)}\right] Fernández and Castro (1985); Herbst and Simon (1978). In the adiabatic regime, the ionization rate by an oscillating monochromatic field is the cycle average of the DC ionization rate corresponding to the instantaneous static electric field at a specific time \mathcal{E}\mathopen{}\mathclose{{}\left(t}\right) Joachain et al. (2011), that is

[TABLE]

As can be seen in the lower panel of Fig. 7, this adiabatic ionization rate corresponds exceptionally well with the average ionization rate obtained from Eq. 21 with $t_{0}=\pi/2\omega$ and $n=1$, i.e., t\in\mathopen{}\mathclose{{}\left[\pi/2\omega;5\pi/2\omega}\right]. It should be noted here that, in this low frequency limit, ECS in both LG and VG, as well as PML in VG, diverge numerically. This phenomenon was briefly discussed above and in more detail in Ref. He et al. (2007). Thus, the low frequency ionization rate has been calculated implementing PML in LG. This can be done without obtaining significant errors for two reasons: (i) we are only interested in the interval $t<5\pi/2\omega$, for which the PML LG method gives a fairly good approximation to the real wave function for low frequencies, and (ii) the ionization rate is not as sensitive to the errors in the wave functions as the error measurement in Eq. 7.

Two further comparisons are made in the lower panel of Fig. 7. The first is an analytical approximation to the low frequency ionization rate. This can be obtained by analyzing the asymptotic behavior of the Airy functions. For low field strengths, \mathopen{}\mathclose{{}\left|\lambda}\right| tends to infinity. In fact, both the real and imaginary part of $\lambda$ tend to $+\infty$ for $\mathcal{E}_{0}\to 0$. Substituting the asymptotic expressions of the Airy functions Abramowitz and Stegun (1972) into Eq. 22, one obtains a polynomial in the electric field multiplied by an exponential function. Solving for the imaginary part of the energy, while retaining only first order terms in the polynomial, then leads to

[TABLE]

The ionization rate of the monochromatic field can then again be obtained by Eq. 23 using Eq. 24 for $\Gamma_{\mathrm{DC}}$. It can be seen to agree with the first two methods for weak fields. The final comparison is made with the expression obtained by Perelomov, Popov, and Terent’ev (PPT) Perelomov et al. (1966) for the adiabatic ionization rate for the zero-range potential in a monochromatic field with amplitude $\mathcal{E}_{0}$

[TABLE]

It again agrees for weaker fields but as opposed to the approximation obtained by the asymptotic analysis it overestimates the ionization rate for strong fields.

VI Long-range Potential

As a final study of the behavior of PMLs in laser-matter interactions, we will look at a long-range potential, namely the one-dimensional ”hydrogen atom”

[TABLE]

This potential does not have the same cut-off spatial behavior as the short-range potential in the previous section, and thus this potential will reach into the absorbing layer. The ground-state energy of the one-dimensional hydrogen atom remarkably comes out as exactly $E_{0}=-1/2\,\mathrm{a.u.}$. First, we study the error as a function of absorption width $d$. This is shown in Fig. 8 for the PML method and in Fig. 9 for ECS. For the PML method, we see similar trends as for the short-range potential with the exception that the errors are much larger. Whereas the PML with $\sigma_{0}=0.001$ and a square absorption function converged to an error of the order $10^{-15}$ for the short-range potential, it converges to around $10^{-6}$ for the long-range potential. It does seem, however, that one is able to obtain lower errors if the absorption width is increased and the absorption function is allowed to increase more slowly. This is indicated in Fig. 8(a), as the nearly singular function will grow slower as $d$ is increased. This leads to an interesting opportunity in implementing a non-uniform grid in the absorption region, significantly reducing the number of grid points required to describe a much wider absorption width. This was done in Ref. Weinmüller et al. (2017) for ECS and was shown to produce excellent results. Figure 10 shows the error as a function of time for the same parameters as for the short-range potential. Here it can be seen that ECS in VG and LG are again indistinguishable at later times. As is evident, the ECS scheme is preferable for a potential that reaches into the absorbing layer.

VII Concluding Remarks

We have presented a simple FDTD scheme implementing PMLs to describe the dynamics of an electron in a laser field. The PML approach has been compared to the well known ECS approach, and we observe that the PML, when implemented in VG, far outperforms the ECS approach when the potential vanishes outside the scaling radius $R_{0}$. Conversely, the traditional ECS approach is much more efficient when implemented for a potential that reaches into the absorption domain. Upon comparing the errors introduced by ECS implemented in the two gauges, little to no difference is observed. On the other hand, when PMLs are implemented in LG, significantly larger errors can be seen. For sufficiently low frequencies, both the ECS LG and VG and the PML VG wave functions blow up, leaving PML LG as the only remaining option of the four methods. Finally, we have demonstrated that the PML implemented in LG is able to reproduce the ionization rate in the adiabatic region, where the other methods considered fail.

Acknowledgements.

Useful comments from Lars Bojer Madsen are gratefully acknowledged. This work was supported by the Villum Kann Rasmussen (VKR) Center of Excellence QUSCOPE. Additionally, H.C.K. and T.G.P. are supported by the Center for Nanostructured Graphene (CNG), which is sponsored by the Danish National Research Foundation, Project No. DNRF103.

APPENDIX A Finite Difference Formulas

The finite difference (FD) approach used in the present paper is based on the Crank-Nicolson scheme Crank and Nicolson (1996). It consists of a combination of the forward (explicit) and backward (implicit) Euler method and reads

[TABLE]

where $t_{j}=j\Delta t$ and H\mathopen{}\mathclose{{}\left(t_{j}}\right) is the Hamilton operator at time $t_{j}$. What remains is to discretize the spatial derivatives. For the PML method, the kinetic term is given by

[TABLE]

Both $c$ and its derivative are known analytically. The derivatives of the wave function are approximated using the second-order approximations

[TABLE]

with \psi_{n}=\psi\mathopen{}\mathclose{{}\left(x_{n}}\right) and $x_{n}=n\Delta x$, where $n$ is an integer running from $-N$ to $N$ for a grid with a total of $2N+1$ equidistant points. These simple FD formulas are one of the advantages of the PML method.

For the ECS method, on the other hand, the transformation leads to modified FD formulas. They can be derived by writing the wave function as a standard approximation using Lagrange interpolating polynomials

[TABLE]

with

[TABLE]

The FD formulas at any particular point are then derived by differentiating Eq. (31) and evaluating the result at said point, all the while keeping in mind Eq. (9). That is,

[TABLE]

where $x_{n}$ is the equidistant grid described above. Here we use $p=1$, which leads to the standard FD formulas for \mathopen{}\mathclose{{}\left|x}\right|<R_{0}. Outside the scaling radius we simply pick up a complex phase factor

[TABLE]

while at the scaling radius we have the non-symmetric formulas

[TABLE]

[TABLE]

As was discussed in Ref. McCurdy et al. (2004), the ECS FD formulas are O\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\Delta x}\right)^{2}}\right] for $x\neq\pm R_{0}$ but only O\mathopen{}\mathclose{{}\left(\Delta x}\right) for $x=\pm R_{0}$. For all calculations in the present paper, we have used $\Delta x=10^{-2}\,\mathrm{a.u.}$ and $\Delta t=10^{-3}\,\mathrm{a.u.}$.