Exterior complex scaling as a perfect absorber in time-dependent problems

Exterior complex scaling as a perfect absorber in time-dependent problems

Armin Scrinzi armin.scrinzi@lmu.de Ludwig Maximilians University, Theresienstrasse 37, 80333 Munich, Germany, EU
and
Wolfgang Pauli Institute, 1090 Vienna
July 13, 2019
Abstract

It is shown that exterior complex scaling provides for complete absorption of outgoing flux in numerical solutions of the time-dependent Schrödinger equation with strong infrared fields. This is demonstrated by computing high harmonic spectra and wave-function overlaps with the exact solution for a one-dimensional model system and by three-dimensional calculations for the H atom and a Ne atom model. We lay out the key ingredients for correct implementation and identify criteria for efficient discretization.

pacs:
42.50.Hz,02.60.Cb,33.20.Xx

I Introduction

The absorption of outgoing parts of the wave function at the boundaries of a finite volume is a key problem for any efficient numerical solution of the time-dependent Schrödinger equation (TDSE) and it has been amply dealt with also in recent literature (see, e.g., Antoine et al. (2008) and references therein). This interest has been renewed in the context of intense laser-matter interactions: speaking in terms of physics, strong fields lead to large ionization and therefore large fluxes out of a central region. For strong field induced electronic and nuclear dynamics in atoms and molecules and high harmonic generation, electrons far from the system play no role and can be disregarded. When solving the TDSE for these processes, one can therefore identify an inner region (a finite volume) where an exact solution is of interest. Outside that region one must, by some means, truncate the solution without compromising the inner region. This is particularly important for higher-dimensional problems involving two or more electrons in order to control the size of the discretization. Out of the large number of approaches towards that goal the majority of computations of strong laser-matter interactions employed one of the following methods: absorbing masks Krause et al. (1992), complex absorbing potentials (CAPs) Riss and Meyer (1996), and exterior complex scaling (ECS)He et al. (2007); Tao et al. (2009).

The two recent numerical studies on ECS have cast doubt on the efficiency Tao et al. (2009) and maybe even the fundamental correctness the method in numerical practice He et al. (2007). In the present paper we will show that ECS indeed is a perfect absorber to full computational accuracy (14 digits). In addition, it allows highly efficient implementation where only a small fraction of the total discretization points are used for absorption. In both respects it far outperforms commonly used monomial CAPs. As a third point, as noted early on McCurdy et al. (1991), ECS is not just an absorber: ideally, it keeps a record of the dynamics in the outer region, which, in principle, could be recovered. We will provide numerical evidence for this fact.

After giving a brief review of ECS, we will present with some care our discretization method, as it plays an important role for correct and efficient implementation of ECS. The general characterization of ECS and a comparison with CAPs is done using a one-dimensional model system, and finally we will present results in three dimensions for the hydrogen atom and a single-electron model of Ne.

Ii TDSE with a laser field

We want to solve the TDSE of the general form

(1)

where will be either a single or three spatial coordinates. then denote and Laplace and Nabla operators, respectively. is a system-dependent binding potential and is the vector potential of the laser field. Here we have chosen the velocity gauge and removed the term by a time-dependent unitary transform. As the initial state we use the lowest energy eigenfunction of the field-free Hamiltonian operator . We use vector potentials with finite duration

(2)

in the time interval with . The peak vector potential is and in 1 and 3 dimensions, respectively. Such pulses with a single or a few oscillations of the electric field, linear polarization, and peak field amplitude at are frequently used as models in numerical studies.

The complete information of the system inside some inner region is contained in the wave function amplitude. For characterizing the accuracy of our results by a single number, we use the overlap between an “exact” solution obtained from a calculation in a very large box and the approximate solution

(3)

The scalar product is restricted to the inner region or a sub-region of the inner region

(4)

and is the corresponding -norm.

A quantity of immediate physical interest is the intensity spectrum of the harmonic response given by the Fourier transform of the “acceleration of the dipole”

(5)

For the comparison, integrals are restricted to the inner region . In general is a highly oscillatory quantity varying by several orders of magnitude. The local error of the spectrum relative to an “exact” spectrum is

(6)

Local averaging in the denominator suppresses spurious spikes due to near-zeros of the spectrum.

Iii Outline of ECS theory

There is a large volume of literature available on complex scaling in general (see, e.g., Balselev and Combes (1971); Aguilar and Combes (1971); Reed and Simon (1982)) and on exterior complex scaling in particular (see, e.g., Nicolaides and Beck (1978); Simon (1979); McCurdy et al. (2004) and references therein). We restrict our summary to communicating the basic idea and to emphasizing the points that are essential for correct numerical implementation. For this we closely follow earlier work, Ref. Scrinzi and Elander (1993). In one dimension, exterior complex scaling consists in continuing the coordinate outside a “scaling radius” into the complex plane

(7)

The effect of the transformation on plane waves at values is

(8)

For positive — outgoing waves to the right side — the functions become exponentially damped with increasing , while for negative they grow exponentially. The corresponding situation with reversed signs arises for . By complex scaling we can distinguish in- from outgoing waves simply by their normalizibility without the need to analyze the asymptotic phase. In a typical discretization we implicitly or explictly use only square-integrable functions, by which we exclude ingoing waves from a complex scaled calculation. This is the key to complex scaling as a perfect absorber: in a well-defined region we have a simple instrument to systematically suppress ingoing waves by just requiring that our solution remain square-integrable. A more elaborate version of this reasoning can be found, e.g., in the appendix of Ref. Scrinzi and Piraux (1998).

The mathematically rigorous theory of complex scaling often uses the alternative point of view that not the wave functions, but rather the operator itself is scaled, while it remains an operator acting on an ordinary Hilbert space of square-integrable functions. We follow this line of reasoning for pointing out a fact of immediate computational relevance. We only give here plausibility arguments and refer the reader to Ref. Scrinzi and Elander (1993) for a more extensive discussion and references to mathematical literature.

One starts from real scaling, i.e. replacing in Eq. (7) with a real number and observes that the transformation

(9)

is unitary. The scaling factor is essential to ensure unitarity. Formally, this transform can just as well be applied to the Hamiltonian by defining

(10)

It is important to note that if is defined on differentiable functions , the transformed operator is defined on the discontinuous functions and its action on these functions is given by

(11)

As a unitary transform leaves the operator’s spectrum unchanged and the scaled dynamics is in a one-to-one relation to the unscaled. Now the hard part of mathematical theory sets in: for a certain class of “dilation analytic” potentials, the operators can be analytically continued to complex values without changes in the bound state spectrum Reed and Simon (1982). The continuous spectrum is rotated around the continuum threshold into the lower complex energy plane by the angle . This is trivial to see for the free particle and the case where the spectrum transforms as

(12)

This property of the continuous spectrum persists when dilation analytic potentials are added and for : the complex scaled Hamiltonian retains a distinct “memory” of the unscaled Hamiltonian. Proof of dilation analyticity can be difficult to find. Beyond some large , where most physical potentials have simple decaying tails, the expected ECS properties are confirmed by numerical experiment.

One can now write an exterior complex scaled TDSE

(13)

Here it is assumed that the potential can be analytically continued to complex values . One may hope that the dynamics generated by (13) is related to the original dynamics, and that for the solution is identical to the unscaled solution . We will demonstrate below that this expectation can be confirmed by numerical experiment to the level of full machine precision.

The main purpose of this brief discussion of ECS theory is to stress the importance of the correct discontinuity in the wave function for defining differential operators. The discontinuity at is intimately linked to the unitarity of the real scaled problem, which in turn secures the conservation of spectral properties of the scaled operator. For given and it has the explicit form

(14)
(15)

The discontinuity condition (15) on the derivative arises from transforming the continuous first derivatives of the original functions. On such functions, one can define the complex scaled Laplacian in analogy to Eq. (11) by “back-scaling” the scaled function , applying the ordinary Laplacian, and forward-scaling the result:

(16)

The factor appears, as the derivative is applied to the back-scaled function rather than to . On continuous functions the scaled Laplacian (and any derivative) is not defined as an operator in the Hilbert space, just as an ordinary Laplacian is not defined on discontinuous functions.

Finally we want to point to the fact that the adjoint operator requires functions with the complex conjugate condition of Eqs. (14,15). This means that for our discretization by a basis set the discontinuities must not be conjugated when going from ket- to bra-vectors. We will show below how this can be easily implemented in a finite element basis.

Iv Discretization

For the discretization of ECS two points are important: (i) the correct implementation of the discontinuity and (ii) good approximation of analyticity. Both can be most conveniently accommodated in a finite element discretization of high rank.

We follow the implementation strategy laid out in Ref. Scrinzi and Elander (1993): each coordinate axis is divided into elements . On each element we choose a set of linearly independent functions that can be transformed to obey the conditions

(17)

We will call the “rank” of the finite element. In principle any set of functions that obeys (17) can be used in a finite-element scheme. In practice we use real-valued polynomials which for enhancing numerical stability we transform to

(18)

with normalization constants . For the element functions (17) Dirichlet boundary conditions are implemented by omitting the functions and . Alternatively, on the leftmost and rightmost intervals we use polynomials times an exponential with and signs on the intervals and , respectively. The conditions on the end element functions are

(19)

The exponent can be adjusted for best performance, but its exact value was found to be uncritical for ECS.

The finite-element ansatz for the total wave function is as usual

(20)

By construction of the , Eqs. (17), continuity across element boundaries is assured by demanding

(21)

Elementwise overlap and Hamiltonian matrices are

(22)
(23)

where denotes the Jacobian function for integration over . The elementwise matrices are added into the overall discretized matrices and such that the last row and column of each elementwise matrix overlaps with the first row and column of the following element (see Fig. 1), which is equivalent to setting the corresponding coefficients equal, Eq. (21). As always in finite element methods, continuity of the first derivative does not need to be imposed (see Ref. Scrinzi and Elander (1993) for a more detailed discussion). and are matrices with

(24)
Figure 1: Placement of the elementwise block in the overall Hamiltonian matrix

For ECS we choose the scaling radii to coincide with the element boundaries . The scaled elementwise Hamiltonian matrices are evaluated by substituting in (23) the Jacobian and the operator with their ECS equivalents

(25)

As we use real functions we could omit the complex conjugation and the resulting matrices are complex symmetric, i.e . The discontinuity (14) is brought into the system by the factor for the integrals : it amounts to multiplying all functions outside the scaling radius by and as the discontinuity does not get complex conjugated, the bra and ket discontinuity factors do not cancel but multiply to . Like in the unscaled case, the continuity condition on the first derivative (15) does not need to be imposed for finite elements. The procedure for constructing the overall matrix is identical to the unscaled case. Replacing by results in the correct (non-hermitian) overlap matrix for the discretized problem. In practice, the matrices and are rarely set up explicitly, as applying the elementwise matrices to the corresponding sections of the coefficient vectors is far more efficient.

There are no specific issues for time propagating the discretized system

(26)

except maybe that very high accuracy was needed for our comparisons. If anything, ECS mitigates the well-known stiffness problems for explicit time-integrators, as high kinetic energies are also associated with large imaginary parts and decay rapidly. We use Runge-Kutta schemes with self-adaptive step size and self-adaptive order up to order 7. Robust error control is achieved by single-to-double-step comparisons. We obtain significant speedups of the propagation by removing states with very high eigenvalues of the field-free Hamiltonian from the simulation space by explicit projection.

V ECS for a 1-d problem

We first investigate ECS for the 1-dimensional “hydrogen atom” with the model potential

(27)

which gives the ground state energy . Here and below use the peak vector potential and the optical period . If interpreted as atomic units, these parameters correspond to peak intensity and wave length . We will show results for FWHM of amplitude of , 5 and 10 optical cycles and total pulse durations of 2,10, and 20 optical cycles, see Eq. (2). The classical quiver amplitude of an electron in this field is atomic units. At the end of a single cycle pulse with this intensity around 20% and after a 5 cycle pulse more than 80% of the electron probability have left the range [-40,40].

Within the framework of this model system we will answer the following questions: Can ECS be considered a perfect absorber? Can the scaling radius be put inside the range of the quiver motion, i.e. ? Does ECS work for long pulses? Which parameters determine the efficiency of ECS? How many discretization coefficients are needed? How does ECS perform compared to monomial CAPs? Does ECS work for length gauge?

v.1 ECS is a perfect absorber

We call an absorber perfect, if the error defined in Eq. (3) is on the level of machine precision. As the “exact” result for comparison we use a unscaled calculation on a large box with a total of discretization coefficients distributed over 120 elements with constant rank . The elements are equidistant except for the infinite end elements and with exponent , Eq. (19).

For ECS we use the parameters and and finite elements that up to are the same as in the unscaled calculation. In the scaled ranges on either end of the axis we use infinite elements and with and exponent . At this point, no attempt was made to minimize the number of coefficients used for absorption by optimizing the scaling parameters. Indeed, with the given parameters we obtain for the errors at the end of the pulses

(28)

The error of the wave function amplitude is about the square root of these values and it remains constant after the initial rise, see Fig. (2). The error level is constant through the whole range and there is a sharp edge to the scaled region, where the wave function is not directly related to the unscaled one. The errors indicate the accuracy limits of our numerical integration scheme and are not determined by ECS. It is therefore fair to say that, at least for the present model, ECS acts as a perfect absorber.

Figure 2: (color online) Evolution of the relative error during a 5-cycle pulse. The denominator is averaged over 5 grid points to avoid spurious spikes. For the pulse parameters and discretization, see text. The sharp rise of error marks the boundaries of the inner region. A plane is drawn at error level (color blue); only a few error peaks in the inner region are above . Away from the center, relative errors are enhanced initially as the wave function is nearly zero.

v.2 Element rank and infinite end elements

The choice of conspicuously high element rank for these very accurate calculations is not by coincidence. Complex scaling depends on analyticity properties of the Hamiltonian. It is therefor not surprising that we observe a strong dependency of the accuracy on the degree to which our discretization can approximate analytic functions. Any localized basis, such as finite elements or B-splines is not analytic by definition because of the finite support of the basis functions. However with increasing polynomial degree, loosely speaking, one gets closer to analytic functions. Table 1 lists the error of the wave function in the range [-35,35] for increasing element rank. As must fall onto element boundaries, we had to choose slightly different values for the different element ranks. The ECS absorption range is discretized with between 36 and 45 exponentially damped functions with such that the sum of coefficients in the scaled and unscaled regions was for all calculations. For the error estimates at each a large box real calculations was performed with same and the same number of points as for ECS in . From Table 1 we see that, depending on the desired accuracy, it is advisable to use polynomial degrees 8 or higher.

=
4 41 40.
5 41 40.
6 41 40.
7 43 39.
9 41 40.
11 41 40.
13 37 42.
15 46 38.
21 41 40.
Table 1: Dependence of the final wave function error on the element rank . All ECS calculations are for a single cycle pulse and a total of discrete coefficients.

For practical purposes we want to mention that the variation of an ECS calculation with and box size is not a safe indicator of its accuracy. We found ECS calculations with fixed , element rank and element sizes to be far better consistent among each other than their error relative to the unscaled result. For reliable accuracy estimates one must vary the scaling radius .

The use of infinite elements at the ends of the simulation box greatly contributes to the good performance of ECS. Table 2 compares a few finite-box calculations with a calculation using infinite end-elements with only 21 discretization points. Only at rather large finite boxes and a larger number discretization points one reaches the infinite element result.

21
20 20
40 40
60 60
80 80
Table 2: Error of ECS calculations with infinite and finite absorption ranges. In all calculations we use a single-cycle pulse, rank and 160 discretization points in . The length of the absorption range is and is the number of coefficients for absorption at each side.

The explanation for this may be as follows: it was noticed in Ref. Tao et al. (2009) that long wave lengths cannot be accommodated in a finite ECS region and deteriorate absorption by reflections. Such long wave lengths have very little structure and should be easily parameterizable. It seems that the exponential tail of our end-element functions is sufficient to accommodate slowly varying long wave-length parts of the ECS wave function. We leave a more detailed investigation of this observation to later work, but conclude that efficient ECS is best done with infinite end intervals.

v.3 Choice of and back-scaling

We find that the quality of the wavefunction in the unscaled region is not affected by the choice of the ECS radius . Table 3 shows the errors for and 40. The general error level in these calculations is slightly higher as we used a lower element rank of in order to be able to make the two elements of the inner region small. The density of discretization points was kept constant through all calculations. We see that the error level is independent of whether the ECS radius is chosen inside or outside the classical quiver amplitude of . Errors only start to rise, when the total size of the box indicated by the number of discretization points becomes too small. This may not be surprising, if we assume that the spatial range of the dynamics remains essentially unchanged by complex scaling: if the box, be it scaled or unscaled, cannot let a particle go the full distance and then return without reflections, distortions must occur.

241 40.
201 20.
160 10.
160 5.
100 10.
80 10.
60 5.
Table 3: Dependence of the final wave function error on ECS radius and on the total number of discretization points .

There is an interesting conclusion that one may draw from this independence of : the fact that significant flux moves into the scaled region and then back out without corrupting the unscaled part of the wave function indicates that also in the scaled region the TDSE dynamics is encoded correctly, although in a different way. Our numerical results are a striking corroboration of this conjecture that was made early on in ECS theory McCurdy et al. (1991). In principle one may hope to recover the unscaled wave function by analytic continuation. This hope for back-scaling, in fact, was the original motivation for introducing the analytic form of functions on the end-elements, as an ordinary finite element function cannot be unambiguously analytically continued. We have not further pursued this possibility for two reasons: first, with larger and larger scaling angles we encountered severe numerical problems, as the back-scaled basis functions become highly oscillatory and cancellation errors destroy the reconstruction of the unscaled wave function; the second reason is the striking success of ECS with just a few points needed for absorption. It is safer and simpler to just discard the small absorption range and use the inner region directly for the evaluation of physical quantities. Yet, if for one reason or another, one wishes to back-scale a time-dependent ECS wave function, our results indicate that such a procedure can be successful. One may in that case use a representation of the scaled region that is less plagued by numerical problems than our exponential basis.

v.4 Choice of scaling angle and exponent

Although with sufficiently large absorption range one can always obtain perfect absorption independent of scaling angle and damping exponent , optimizing these parameters in a given situation allows to obtain good absorption with very few absorption points. Figure 3 shows the error for n=1 and n=5 cycle calculations with and 20 absorption points on either end of the interval. The exact choice of the parameters is not critical for the calculations, where full accuracy is reached for rather large parameter ranges. As to be expected, the 5-cycle calculation with is most sensitive to and , but still in a range of and around the optimal values accuracy deteriorates only by 2 orders of magnitude to the still acceptable value of . There is a clear anti-correlation between and , which may be explained looking at the oscillatory behavior of the back-scaled exponential . We conjecture that the effective back-scaled wave number is the relevant parameter for efficient absorption. Correlation between the parameter and nearly vanishes and optimization can safely be performed for each parameter independently.

Figure 3: (color online) Error as a function of scaling angle and exponent for n=1 and n=5-cycle pulses. For the 5 cycle pulse, longer wave length reach the boundaries; optimal exponent and scaling angle are smaller and a longer (20 point) absorption range significantly increases accuracy.

v.5 Comparison to complex absorbing potentials

A popular and comparatively straight forward way of absorbing outgoing flux are complex absorbing potentials (CAPs). The basic idea is to add at the end of the simulation box a potential with a negative imaginary part, which leads to exponential damping of the wave function. In this simplest form, the method can be considered a differential form of absorption by mask functions, where at preset intervals a certain part of the wave function is removed. A fundamental limitation of these methods is that they — even in principle — cannot be strictly reflectionless. The attempt to obtain minimal reflections has lead to range of models, partially including real parts into the potential and adjusting to specific physical situations (see, e.g. Muga et al. (2004)).

It is beyond the scope of the present work to perform a comprehensive study of CAP for the present type of problems. Rather, we use the simple and well-investigated monomial CAPs Riss and Meyer (1996)

(29)

for polynomial degrees with optimized in each calculation. The criterion for our comparison with ECS is the number of discretization points needed for a given level of absorption.

Results are shown in Table 4. With a finite absorption range, ECS outperforms CAP roughly by one or two orders of magnitude. However, when we use infinite end elements with ECS (discretized by only 21 points), we can reach absorption to machine precision. We could not find a similar adjustment for CAP.

Method or q
ECS 21 0.6
ECS 20 10 0.6
ECS 40 20 0.5
CAP 20 10 4
CAP 20 10 6
CAP 40 20 4
CAP 60 30 4
Table 4: Accuracy of ECS and CAP for different absorption ranges and number of absorption coefficients . Scaling angle and absorption strength for ECS and CAP, respectively, were optimized. The errors are calculated at the end of a single-cycle pulse.

v.6 High harmonic spectra

Although the error is a meaningful measure for wave function accuracy, it cannot be directly related to the error of a given observable. Figure 4 shows the accuracy of ECS high harmonic spectra of 1- and 5-cycle pulses relative to a real calculation. We find errors on the level between and and we could not get much better agreement than this irrespective of discretization and scaling parameters. Again this error is related to the numerical limits of our discretization and propagation schemes: the wave function error is and the high frequency signal is suppressed by relative to the fundamental peak, making a relative error of the suppressed signal of the order quite plausible. Indeed we find similar errors when comparing different, but equally accurate purely real calculations. More disquieting is the error at the fundamental frequency, which does not appear in large box real calculations. We were not able to locate the origin of this error: it persists through variations of , specific discretizations, different time-discretizations, and also for the 3-d hydrogen calculation below (cf. Fig. 5). The error appears to be related to an artificial overall modulation of the signal by the driving field, possibly related to internal normalizations during propagation. Note that by construction normalization errors do not appear in the wave function accuracy measure , Eq. (3). We believe, however, that this error is acceptable for all practical purposes.

Figure 4: (color online) High harmonic power spectrum for a 5-cycle pulse (upper panel). Lower panel: Accuracy of the high harmonic spectrum with different ECS parameters and discretizations. Curve A is the error for relative to a fully converged real calculation. The choice of has has the largest influence on accuracy: curve B, the difference between two calculations with and closely follows the overall error curve A. At fixed the influence of the other ECS parameters and discretization is small: curve C compares a calculations using =(40,20,0.7,0.7) and finite element rank with (40,40,0.5,0.3) and rank 41.

v.7 ECS fails in length gauge

For field-interaction in length gauge

(30)

ECS completely fails in the time-dependent case. The reason for this behavior was pointed out in Ref. McCurdy et al. (1991): when length gauge Volkov solutions are complex scaled their asymptotic behavior becomes dependent on the sign of the field strength and alternates between damping and growth. The convenient distinction between in- and outgoing waves by their norms is lost. In the language of mathematical theory, is not a dilation analytic potential, and severely so: complex scaling transforms the spectrum of the Stark problem from purely continuous into purely discrete and all discrete eigenvalues of the scaled Stark Hamiltonian have imaginary parts Herbst (1980). This is sharp contrast to dilation analytic potentials where the bound state energies remain unchanged and the continuous spectrum is only rotated into the lower complex plane.

Vi Calculations for H and model Ne

In order to demonstrate the applicability of ECS to realistic problems, we show calculations for the H atom with

(31)

and a single electron model of the Ne atom with the potential

(32)

We use the parameters given in Table 5, for which our model reproduces the ground and first few excited state energies of Ne.

1 -1 0
2 -0.3 0.5
3 -2.05 2
4 1.23 1
Table 5: Parameters for the Ne model potential Eq. (32)

We use linearly polarized pulses with the same pulse shape and peak intensity as in the preceding section and pulse durations of 1 and 10 optical cycles. The calculations are done in polar coordinates with a spherical harmonics basis on the angular coordinates and high rank finite elements on the -coordinate. Again an infinite last element is used.

There are no surprises: convergence patterns and accuracy are very similar to the one-dimensional model. Figure 5 shows the harmonic spectrum for hydrogen at 1 cycle together with errors for different ECS and discretization parameters. Error estimate here is by comparison to an calculation.

Figure 5: (color online) High harmonic power spectrum from the hydrogen atom for a 1-cycle pulse (upper panel). Lower panel: Error with ECS parameters relative to a calculation (curve A). The relative difference to a calculation with (40,40,0.4,0.4), curve B, underestimates the error. The calculation is converged with 20 angular momenta on the given level of accuracy. More angular momenta do not change the result. At 15 angular momenta (curve C) accuracy deteriorates.

No new problems appear due to the more general Ne model potential (32). Figure 6 shows high harmonic spectra from a H and Ne for a 10-cycle pulse. Accuracy estimates were obtained by varying ECS radius .

Figure 6: (color online) High harmonic power spectra from the hydrogen and a neon model with a 10-cycle pulse (upper panel). Relative accuracies shown in the lower panel are somewhat poorer with the longer pulse, in particular for Ne where due to the higher ionization potential the signal is very weak.

Vii Discussion

As we find high numerical stability and excellent performance of ECS as an absorber, the questions arises what are the reasons for the numerical problems reported in Refs. He et al. (2007); Tao et al. (2009), where ECS was applied to essentially the same systems. One obvious source of inaccuracies may lie in possibly low order discretizations. Unfortunately, in neither publication an investigation of the dependence of the observed effects on discretization is shown.

Certainly the choice of finite box-sizes lowers the performance, but according to Table 2 with the very large absorption ranges of 80 Bohr used in Ref. Tao et al. (2009), excellent results should be achievable.

A possible source of the observed difficulties may be the treatment of the overlap matrix. As discussed above we replace the ordinary overlap by the pseudo-overlap matrix . With this choice and as we use strictly real finite element functions we obtain complex symmetric matrices for zero field and . There are no explicit statements about in Refs. He et al. (2007); Tao et al. (2009). Usually, finite difference methods imply (an approximation to) the identity operator for overlap. The B-spline method used in Tao et al. (2009) requires a choice for and Eq. (20) of Ref. Tao et al. (2009) seems to imply that indeed the identity was used as an overlap matrix.

The comment on the non-orthogonality of the eigenvectors of the non-normal scaled Hamiltonian in Tao et al. (2009) also seems to indicate, that the ordinary, unscaled overlap matrix was used. Clearly, the eigenvectors of the eigenproblem

(33)

will not be orthogonal in general. However, we find that all eigenvectors of the complex scaled generalized eigenproblem

(34)

are pseudo-orthogonal and can be normalized in the sense

(35)

Then the matrix has a diagonal representation

(36)

and the spectral values appear as discrete approximations to the true ECS spectrum with strictly non-positive imaginary parts . We do not have mathematical proof for this property of the discrete complex scaled system, but we find it valid in all our calculations on the level of computational accuracy. If on the other hand we use the ordinary overlap matrix , we invariably obtain a few eigenvalues with which will cause long-term instability of the time-propagation. Possibly, this is the ultimate reason for the numerical instabilities observed in Refs. He et al. (2007); Tao et al. (2009).

Viii Summary

We have demonstrated that ECS can serve as a perfect absorber of outgoing flux in the sense that in the unscaled inner region it exactly matches a purely real calculation on a sufficiently large grid. We were able to push the agreement to relative error of . The corresponding errors in the wave function amplitude are . Both errors are at the limits of our numerical integration scheme.

Furthermore we have evidence that also in numerical practice ECS does not just act as an absorber but conserves dynamical information during excursions into the absorbing region: even when the quiver motion takes flux deeply into the “absorbing” region, the returning flux is identical to the flux in a purely real calculation.

For this, we found the following points essential:

  • implementation of the correct scaled derivatives, including bra-functions with unconjugated discontinuity,

  • the use of “infinite” absorption ranges, which we discretized by polynomials times an exponential,

  • the use of high rank discretization also in the inner region to reach the highest accuracies.

Point (i) leads to a complex symmetric, in particular not positive definite discrete overlap matrix which must not be approximated by a positive definite matrix.

Following these rules, we encountered no numerical difficulties in the inner region or in the absorbing region, using a standard explicit Runge-Kutta scheme for time integration. As a tendency, large scaling angles favor good absorption, in many cases we used , which corresponds to an almost purely imaginary continuous energy spectrum . In our basis we found the scaling angle ultimately to be limited by numerical instabilities due to the complex symmetric overlap matrix. As excellent absorption can be achieved with as few as 20 discretization coefficients in the absorbing region, pushing the scaling angle to the numerical limits is not necessary in general and scaling angles of served well for our purposes. In general, we found the scheme numerically robust and not very sensitive to the scaling parameters. The option of back-scaling the solution to was abandoned due to severe cancellation errors in the related transformations.

When judging the accuracy of an ECS calculation, it is important to vary the ECS radius . Our comparison with a real calculation indicates the variation of the result with different gives realistic error estimates. Other parameters such as rank of the finite elements, length of the absorption range, or scaling angle are of secondary importance.

ECS in the present implementation outperforms simple monomial CAPs. ECS errors were one or two order of magnitude smaller than CAP errors with the same spatial discretization. When using infinite end intervals, the advantage of ECS can even reach 12 orders of magnitude! We are aware of the fact CAPs can be greatly improved by a variety of measures (see, e.g., Muga et al. (2004)). However, in general these require tuning of the CAP parameters to a given situation. Even with that, we do not expect to reach comparable accuracies with CAPs as we could demonstrate for ECS.

We believe that ECS solves the absorption problem for the present class of systems in any discretization, where implementation of (i)-(iii) is possible. Recovery of asymptotic information, such as electron spectra, may be inherently difficult as the total amount of information that is contained in the scaled discretization is too small, which manifests itself in the ill-conditioning of the back-scaling problem. However, in Ref. Caillat et al. (2005) we presented a scheme for computing electron spectra under the assumption of a perfect absorber, at that point formulated for CAPs. An adaptation of that scheme to ECS and extension to few-body dynamics will be investigated in future work.

Ix Acknowledgement

This work was supported by the Viennese Science Foundation (WWTF) via the project ”TDDFT” (MA-45).

References

  • Antoine et al. (2008) X. Antoine, A. Arnold, C. Besse, and A. Ehrhardt, Matthiasand Schädle, Comm. Comp. Phys. 4, 729 (2008).
  • Krause et al. (1992) J. L. Krause, K. J. Schafer, and K. C. Kulander, Phys. Rev. A 45, 4998 (1992).
  • Riss and Meyer (1996) U. V. Riss and H.-D. Meyer, J. Chem. Phys. 105, 1409 (1996).
  • He et al. (2007) F. He, C. Ruiz, and A. Becker, Phys. Rev. A 75, 053407 (2007).
  • Tao et al. (2009) L. Tao, W. Vanroose, B. Reps, T. N. Rescigno, and C. W. McCurdy, Phys. Rev. A 80, 063419 (2009).
  • McCurdy et al. (1991) C. W. McCurdy, C. K. Stroud, and M. K. Wisinski, Phys. Rev. A 43, 5980 (1991).
  • Balselev and Combes (1971) E. Balselev and J. Combes, Comm. Math. Phys. 22, 280 (1971).
  • Aguilar and Combes (1971) J. Aguilar and J. Combes, Comm. Math. Phys. 22, 269 (1971).
  • Reed and Simon (1982) M. Reed and B. Simon, Methods of Modern Mathematical Physics (Academic, New York, 1982), p. 183ff.
  • Nicolaides and Beck (1978) C. Nicolaides and D. Beck, Phys. Lett. A 65, 11 (1978).
  • Simon (1979) B. Simon, Phys. Lett. A 71, 211 (1979).
  • McCurdy et al. (2004) C. McCurdy, M. Baertschy, and T. Rescigno, J. Phys. B 37, R137 (2004).
  • Scrinzi and Elander (1993) A. Scrinzi and N. Elander, J. Chem. Phys. 98, 3866 (1993).
  • Scrinzi and Piraux (1998) A. Scrinzi and B. Piraux, Phys. Rev. A 58, 1310 (1998).
  • Herbst (1980) I. W. Herbst, Comm. Math. Phys. 75, 297 (1980).
  • Muga et al. (2004) J. Muga, J. Palao, B. Navarro, and I. Egusquiza, Physics Reports 395, 357 (2004).
  • Caillat et al. (2005) J. Caillat, J. Zanghellini, M. Kitzler, O. Koch, W. Kreuzer, and A. Scrinzi, Phys. Rev. A 71, 012712 (2005).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
176897
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description