An optimized absorbing potential for ultrafast, strong-field problems

To theoretically describe highly nonperturbative interactions — such as strong-field physics — in a fully quantitative manner, the best option is usually to numerically solve the time-dependent Schrödinger equation (TDSE). One of the most popular approaches to practically solving the TDSE represents the wave function on a finite spatial grid with boundary conditions applied at its edges. In general, such a grid needs to be large enough so that there are no reflections from the boundaries which behave as infinitely hard walls. Otherwise, the reflections might lead to unphysical changes in the observables. For example, an ionized wavepacket reflected from the boundary back to the origin might be driven by the laser field into bound states, thus reducing the total ionization yield.

The most direct way to avoid such spurious reflections is to move the boundary further away. Since the grid density is fixed physically by the highest energy, however, this requires more grid points which, in turn, incurs a greater computational cost. In fact, the large grids required to describe current experiments have become a key bottleneck to improving the numerical efficiency of solving the TDSE, especially as laser wavelengths push beyond 800 nm.

Fortunately, if the wave function at large distances can easily be reconstructed or is not of interest, it can be absorbed at a sufficiently large distance that it does not affect the dynamics at small distances. Applying such absorption techniques, one can generally reduce the box size significantly. The absorb-and-reconstruct strategy was probably first developed by Heather and Metiu Heather and Metiu (1987) which they demonstrated for strong-field dissociation. Their work has been adopted in hundreds of papers since. A new implementation following this philosophy Tao and Scrinzi (2012); Morales et al. (2016) has proven similarly effective.

Among the various methods to effect absorption at the boundary, the most widely used — and probably the simplest — method is the complex absorbing potential (CAP) Kosloff and Kosloff (1986); Vibok and Balint-Kurti (1992a, b); Macias et al. (1994); Riss and Meyer (1996); Ge and Zhang (1998); Hussain and Roberts (2000); Manolopoulos (2002); Poirier and Carrington (2003); Muga et al. (2004); Tarana and Greene (2012); Strelkov et al. (2012); Yue and Madsen (2013); Krause et al. (2014) or the closely related masking function Krause et al. (1992). Another increasingly popular absorbing-boundary technique is exterior complex scaling Rescigno et al. (1999); Rescigno and McCurdy (2000); He et al. (2007); Tao et al. (2009); Scrinzi (2010), where one rotates the coordinate into the complex plane at large distances. Other, less common, methods to treat the boundary reflection include time-dependent coordinate scaling Sidky and Esry (2000); Hamido et al. (2011); Miyagi and Madsen (2016), the interaction representation Zhang (1990); Williams et al. (1991); Yao and Chu (1992), and Siegert-state expansions Yoshida et al. (1999). While these methods are local in time and vary from exact to approximate, it is also possible to construct a perfectly transparent boundary using Feshbach projection techniques Muga et al. (2004). The disadvantage of such methods is that they require the wave function from previous times and are thus nonlocal in time. In this work, we will focus on the CAP due to its popularity and the simplicity of its implementation. Our goals are to make it both more efficient and more effective.

Although the CAPs utilized in previous studies are predominantly polynomials Riss and Meyer (1996); Ge and Zhang (1998); Poirier and Carrington (2003); Muga et al. (2004); Tarana and Greene (2012); Yue and Madsen (2013); Krause et al. (2014), other types of absorbing potentials such as $\cos^{2}$ Strelkov et al. (2012), Pöschl-Teller ($1/\cosh^{2}$) Kosloff and Kosloff (1986), and a pseudo-exponential [$\exp(-x^{-n})$] Vibok and Balint-Kurti (1992a, b) have also been used. In most of these papers, the CAP’s performance is examined by studying the dependence of the reflection $R$ and transmission $T$ coefficients on the energy. Riss and Meyer Riss and Meyer (1996), for instance, carefully investigated the properties of $R$ and $T$ for polynomial CAPs, finding some difficulty in treating low energies. They characterized their optimized potential parameters in terms of the absorbed energy ratio $E_{\rm max}/E_{\rm min}$, where $E_{\rm min}$ and $E_{\rm max}$ indicate the minimum and maximum energies, respectively, for which absorption exceeds a given value. The maximum ratio they considered, 30, is too small, however, for typical strong-field electronic dynamics. We will, for instance, consider $E_{\rm max}/E_{\rm min}$=500. Vibok and Balint-Kurti Vibok and Balint-Kurti (1992a, b) presented a more optimal CAP — the $\exp(-x^{-n})$ type — for heavy particles, but the range of absorbed energies was still insufficient for strong-field problems.

Even though $R$ and $T$ provide a clear, quantitative measure of performance, studies of CAPs in strong-field problems utilizing them can hardly be found. Their absence is likely due to the inherent time-dependent nature of the strong-field problem and the authors’ consequent focus on wavepacket behavior, losing track of the connection with $R$ and $T$. In contrast, we will adopt $R$ as the figure of merit for designing our absorbing potentials for the strong-field problem, incorporating it as a critical piece in our systematic CAP construction method.

The major advantage of the CAP is its simplicity. The major disadvantage is that it has required a relatively large spatial range to be effective, thus consuming non-negligible computational resources. In this paper, we improve the performance of the CAP and systematically design a more optimal — yet general — CAP for strong-field processes. Specifically, we provide an optimized CAP with a factor of 3–4 reduction in the absorption range compared to some widely-used CAPs Muga et al. (2004). Our optimized CAP absorbs at a prescribed level over a large enough energy range to be useful for strong-field processes.

To be clear, while we optimize our CAP for the strong-field problem, it can be readily adapted and re-optimized for other problems following the procedures we outline below.

Since a time-dependent wavepacket can always be written as a superposition of the time-independent scattering states, we use the time-independent reflection coefficient as a quantitative tool for characterizing and designing an optimal CAP. We will require the CAP’s reflection coefficient to remain below a cutoff value $R_{c}$, $R\leqslant R_{c}$, over a given energy range $E_{\rm min}\leqslant E\leqslant E_{\rm max}$. Since the spatial region devoted to the CAP near the edge of the grid is unphysical, we wish to minimize the computing resources it consumes as much as possible. Therefore, in this work, our optimization efforts focus on reducing the absorption range $x_{R}$, as defined in Fig. 1, while meeting the absorption criteria above.

We study one-dimensional CAPs since they can be easily adapted to higher dimensions, obtaining the reflection coefficient $R$ by solving the Schrödinger equation,

as indicated schematically in Fig. 1. We consider the potential $V(x)$ to be one of the CAPs listed in Table 1. The shapes of all the CAPs considered are generically as in Fig. 1 and are controlled by the following parameters: $\alpha^{2}$ is the strength of the potential, $\beta$ mainly determines its width, and $x_{0}$ is a shift. These are the parameters that will be varied to optimize the CAPs.

The JWKB-based CAP obtained by Manolopoulos Manolopoulos (2002) — labeled M-JWKB in Table 1 — is qualitatively different from the others, however, in that it requires no optimization. This simplicity is certainly one of its strengths and derives from the fact that its reflection coefficient effectively decreases monotonically from unity at zero energy to $e^{-\sqrt{2}\pi/\delta}$ at infinite energy. Its explicit expression is in terms of the Jacobi elliptic function $\text{cn}(u,k)$,

with $x_{0}$=$2.622/(2\delta k_{\text{min}})$ where $k_{\text{min}}$=$\sqrt{2mE_{\text{min}}/\hbar^{2}}$ Manolopoulos (2002). One simply chooses $\delta$ from the condition $R(E_{\text{min}})=R_{c}$.

The first five CAPs in the table are defined to be non-zero only for $0\leqslant x\leqslant x_{0}$ and to vanish identically for $x>x_{0}$. The remaining CAPs are defined for all $x$, but vanish exponentially with $x$. The first four CAPs are some of the most commonly used, with the cosine masking function recast as a CAP using $e^{-iV(x)\Delta t}\sim\cos^{1/8}[(x-x_{0})/\beta]$.

We include the Pöschl-Teller potential because it is well known to be reflectionless for specific real values of $i\alpha^{2}$, suggesting that it might have advantageous properties as a CAP. It can be shown analytically, however, that this property no longer holds for complex $i\alpha^{2}$. In the process of optimizing it for the present purposes, we found that shifting its center off the grid minimized $x_{R}$, leaving only its exponential tail on the grid. This result suggested using instead the simpler family of exponential CAPs included in the table.

We calculate the reflection coefficient numerically using the finite-element discrete-variable representation (FEDVR) McCurdy et al. (2001); Schneider et al. (2006) and eigenchannel R-matrix method Aymar et al. (1996). The reflection coefficient can also be calculated analytically for several of the potentials in Table 1. However, we give the analytic solutions (derivations in the appendices) only for the CAPs we propose — namely, the single-exponential and double-exponential CAPs. The double-sinh potential has no analytic solution to the best of our knowledge. In these cases, we confirmed that the R-matrix reflection coefficients agreed with the analytical $R$ to several significant digits.

Since our primary goal is to systematically design an absorbing potential for the strong-field ionization problem with predetermined properties, we will use atomic units for the remainder of our discussion. Our results can be readily applied to other problems, though, using the derivations in the appendices in which the masses and SI units are explicitly retained.

We demonstrate our optimization procedure in detail below first for the single-exponential CAP since it is the simplest to optimize. It also establishes a few key results important for the optimization of our recommended CAPs: the double-exponential and double-sinh potentials. Whether the solution is analytical or numerical, the procedure we describe for optimization is the same and can be applied to other CAPs as well. In fact, this is what we have done for the comparison in Sec. IV.

The values of $R_{c}$, $E_{\text{min}}$, and $E_{\text{max}}$ that we will focus on for this discussion are

We chose this energy range to cover $0.1\hbar\omega\leqslant E\leqslant 14U_{p}$ for an 800-nm laser pulse at $10^{14}\,\text{W/cm}^{2}$ ($U_{p}$ is the pondermotive energy: $U_{p}=I/4\omega^{2}$ with $I$ the intensity and $\omega$ the laser frequency). This energy range includes essentially all photoelectrons one would expect to be produced in this typical pulse. In momentum, which is more convenient for the analytical $R$, this range corresponds to

Note that 14$U_{p}$ exceeds the highest-energy electrons one would normally expect in a strong-field problem, but we will show below that this choice has little effect on the resulting $x_{R}$. Finally, we use $V_{c}$=$10^{-4}\,\text{a.u.}$ to define $x_{R}$ from $|V(x_{R})|$=$V_{c}$.

To determine which CAP — among those listed in Table 1 — is the best, we numerically searched for their optimal parameters, assuming they are purely imaginary potentials. From the discussion above, we know that each could be improved by including a real potential and a complex boundary condition, but we expect — and confirmed with spot tests — that the relative performance of the CAPs will remain qualitatively the same. As mentioned previously, we selected the CAPs to compare based on their apparent popularity in the literature or on the claims made for their performance.

In optimizing these CAPs, we follow the principles described in previous sections that the width of the CAP determines the long-wavelength absorption; and the depth, the short-wavelength. The optimization is then reasonably straightforward once we identify the parameters corresponding to the width and depth.

In Table 2, we list the optimal parameters we have found for our $E_{\text{min}}$, $E_{\text{max}}$, and $R_{c}$. The table includes the resulting values of $x_{R}$, and we expect that they are the optimal values to within a few percent. Note that we used $\delta$=0.1 for the M-JWKB CAP based on the solution of $R(E_{\text{min}})=R_{c}$ taken from Fig. 3 of Ref. Manolopoulos (2002). We show in Fig. 10 the corresponding reflection coefficients.

The cosine masking function should be regarded as a polynomial CAP since its behavior in $0\leqslant x\leqslant x_{0}$ for the optimal parameters of Table 2 is largely indistinguishable from the quadratic CAP — thus its $x_{R}$ is identical to the quadratic CAP. Similarly, for the optimal parameters we found, only the exponential tail of the Pöschl-Teller potential remained on the grid, making its performance essentially identical to that of the purely imaginary single-exponential CAP.

The absorption ranges $x_{R}$ listed in Table 2 display a surprisingly large range — more than a factor of 4. Comparing only the purely imaginary potentials, the exponential and Pöschl-Teller forms are more than 30% more efficient than the quadratic CAP. They are also more efficient than the quartic CAP. So, while the exponential form generally seems better, the majority of the disparity in Table 2 arises from adding a real part to the CAP and imposing a complex boundary condition.

From Table 2, the best performance is found for the double-exponential and double-sinh CAPs, outperforming the next-best CAP — the complex-valued, complex-boundary-condition, single-exponential CAP — by roughly 40%. Compared to the next-best purely imaginary CAP, they hold nearly a factor of 3 advantage in $x_{R}$. For reference, we tested the strategy of adding a real part and a complex boundary condition to the quadratic CAP and found $x_{R}$ shrank only to about 70 a.u. So, while pursuing this strategy with the other CAPs in the table would reduce their $x_{R}$, we believe the double-exponential and double-sinh CAPs would remain the best. Interestingly, since the de Broglie wavelength at $k_{\text{min}}$ is 57 a.u., our best CAP manages its efficient absorption in a range of only about half this longest wavelength.

Our recommendation, therefore, is to use the double-sinh CAP when its singularity at the boundary causes no numerical difficulties. In the cases that it does, then the double-exponential CAP is the best choice. The remainder of our discussion will thus focus on these two CAPs.

The discussion and optimization so far has centered on the $E_{\text{min}}$, $E_{\text{max}}$, and $R_{c}$ from Eq. (3). Other choices may well be needed, however, for other choices of laser parameters or calculational goals. We thus present in this section the optimal parameters for the double-exponential and double-sinh CAPs for a selection of likely changes in $E_{\text{min}}$, $E_{\text{max}}$, and $R_{c}$.

To verify that the improved performance of our recommended CAPs does indeed carry over to the time-dependent problem and its numerical solution, we solve the TDSE for free-electron wavepacket propagation. The wavepacket we use possesses a broad momentum distribution comparable to the target momentum range from Eq. (4), as shown in Fig. 12(a).

We again use FEDVR as the spatial representation and propagate the wave function using the short-time evolution operator

where the Hamiltonian includes the CAP $V(x)$,

and $H_{0}$ is merely the kinetic energy in the present case. We evaluate $e^{-iH\delta}$ using the split-operator form Schneider et al. (2011)

Since $V$ is diagonal in FEDVR, $e^{-iV\delta/2}$ can be easily evaluated and applied. Moreover, in this form, the singularity in the double-sinh CAP causes no problems whatsoever. The action of the remaining term in $H_{0}$ is calculated via a Padé approximation Blanes et al. (2009).

Equation (29) is a simple and convenient way to add a CAP to any propagator. In many cases, the alternative, keeping the CAP in $H$, would require modifications of the propagation algorithm or parameters to handle its non-Hermiticity or the singularity of the double-sinh CAP — or both. These issues were discussed somewhat in Sec. III.3 and more in Ref. Manolopoulos (2002). Using Eq. (29) avoids these concerns and is more than sufficient for the application of the CAP.

The FEDVR element distribution is uniform in the range $-600\,{\rm a.u.}\leqslant x\leqslant 1200\,{\rm a.u.}$, and we require $\psi=0$ at the boundaries. Given that the wavepacket is initially centered near $x=0$, this box is large enough for the wavepacket to propagate for 10 fs without significant reflection at the boundaries even without CAPs. This will be our reference solution. We carry out a second, identical propagation but apply the double-sinh CAP at $536\,\text{a.u.}\lesssim x\leqslant 600\,\text{a.u.}$. For this example, we choose the CAP designed to have $R_{c}=10^{-6}$ using the parameters shown in Table 4. We thus compare the wavepacket with and without applying the CAP. All the results have been tested to be converged to at least 3 digits with respect to all numerical parameters.

Figure 12 shows the two wavepackets at different times. It is clear that the CAP is performing as expected since the wavepacket decays in the absorbing region without any of the characteristic oscillations of reflections visible — at least without enlarging the plot by a factor of four or five. For comparison, the wavepacket without a CAP equally clearly shows the reflection oscillations at $t$=13.1 fs for reflections from the boundary at $x$=1200 a.u.

In addition, the with-CAP wavepacket is numerically unaffected before entering the absorbing region, agreeing with the without-CAP wavepacket to at least 3 digits for $x\lesssim 536$ a.u. for all times, even as more than 70% of the wavepacket has been absorbed. This agreement shows that the absorption range $x_{R}$ defined in the time-independent study is consistent with the results from the time-dependent propagation.

Enlarging Fig. 12 by a factor of at least four or five will reveal the tiny reflection ripples in the with-CAP wavepacket near and in the absorption region. Their relative magnitude is about $10^{-3}$=$\sqrt{R_{c}}$ as expected. Per the discussion in Sec. V.2, such oscillations are unavoidable with a CAP and testing with other CAPs and values of $R_{c}$ further support the conclusions there. For example, the oscillation for $R_{c}$=$10^{-3}$ becomes fairly noticeable, which suggests that $R_{c}$ should in practice be no larger than $10^{-4}$ — i.e., two digits in the wave function — to provide quantitatively reliable results. Finally, we note that the roughly 15-fs propagation time is comparable to a typical strong-field calculation, bringing some realism to this simple demonstration.

We have presented a systematic study to boost the performance of complex absorbing potentials. Based on the time-independent reflection coefficient, we were able to quantitatively design the most optimal absorbing potential of a given form. In particular, for ultrafast, strong-field TDSE solvers, the optimal CAP parameters should be determined by the absorbing energy range required for the laser parameters and by the desired accuracy of the TDSE solutions.

We proposed two new CAPs — namely, the double-sinh CAP and the double-exponential CAP — that significantly outperform the CAPs currently in standard use. Their superiority was demonstrated through comparison with optimized versions of most of the CAPs commonly found in the literature. Both of our proposed CAPs overcome the primary impediment to efficient performance — absorption of long wavelengths — while also absorbing a large range of energies that covers basically all strong-field processes. Because we quantified the CAP’s performance and identified $x_{R}$ as the figure of merit for their efficiency, we could show that using an exponential CAP already improved on the common quadratic CAP’s performance by one third. A further factor of almost three was gained, however, by adding a real part to the CAP — a strategy well known in other uses of CAPs, but not in strong-field applications.

We highly recommend the double-sinh CAP for local spatial representations, such as FEDVR where the potential is diagonal. It is efficient, easy to use, and easy to adapt to different absorption criteria — i.e., energy range and level of absorption. Incorporating it into the time propagation via split-operator methods is easy and effective.

For other spatial representations, the double-sinh and the double-exponential CAPs are equally recommended. However, care might need to be taken for the double-sinh singularity close to the boundary. In case the double-sinh singularity is a problem for the propagator, the double-exponential CAP should be chosen. Although optimization of the double-exponential CAP is more involved than for the double-sinh CAP, it is still fairly straightforward. Its optimization procedure, along with that for the double-sinh CAP, is detailed in this work.