An efficient multigrid calculation of the Far Field map for Helmholtz and Schrödinger equations
Abstract
In this paper we present a new highly efficient calculation method for the far field amplitude pattern that arises from scattering problems governed by the dimensional Helmholtz equation and, by extension, Schrödinger’s equation. The new technique is based upon a reformulation of the classical realvalued Green’s function integral for the far field amplitude to an equivalent integral over a complex domain. It is shown that the scattered wave, which is essential for the calculation of the far field integral, can be computed very efficiently along this complex contour (or manifold, in multiple dimensions). Using the iterative multigrid method as a solver for the discretized damped scattered wave system, the proposed approach results in a fast and scalable calculation method for the far field map. The complex contour method is successfully validated on Helmholtz and Schrödinger model problems in two and three spatial dimensions, and multigrid convergence results are provided to substantiate the wavenumber scalability and overall performance of the method.
1Introduction
Scattering problems are of key importance in many areas of science and engineering since they carry information about an object of interest over large distances, remote from the given target. Consequently, ever since their original statement a variety of applications of scattering problems have arisen in many different scientific subdomains. In chemistry and quantum physics, for example, virtually all knowledge about the inner workings of a molecule has been obtained through scattering experiments [31]. Similarly, in many reallife electromagnetic or acoustic scattering problems information about a far away object is obtained through radar or sonar [13], intrinsically requiring the solution of 2D or 3D wave equations.
The near and far field of a scattering problem are specific regions in space identified by the distance from the scattering object. In these regions the scattering solution is subject to certain physical and mathematical properties. In the far field, the solution behaves as a spherical outgoing wave with an amplitude that decays as , where is the distance to the scattering object. In this regime, the (asymptotic) solution can be written as a product of an angular part and a radial part. The near field is the region just outside the scattering object, where the solution has not yet reached its asymptotic form. In electromagnetic scattering, for example, the near field is dominated by dipole terms. Predicting the near and far field solution is important for many presentday industrial applications (e.g. radar), but also plays a key role in physical measurement systems such as tomography, near field microscopy and MRI.
However, far field maps are not limited to 2D or 3D scattering problems. New stateoftheart experimental techniques in physics measure the full 4 angular dependency of multiple fragments escaping from a molecular reaction [30]. Through these experiments, the reaction rates involving multiple fragments can be detected in coincidence. Many experiments are being planned at facilities such as e.g. the DESY Freeelectron laser (FLASH) in Hamburg or the Linac Coherent Light Source (LCLS) in Stanford. The accurate prediction of the corresponding amplitudes starting from first principles requires the use of efficient numerical methods to solve the highdimensional scattering problems, which can scale up to 6D or 9D in this context. Indeed, after discretization one generally obtains a large, sparse and indefinite system of equations in the unknown scattered wave. Direct solution of this system is usually prohibited due to the massive size of the problem in higher spatial dimensions. Preconditioned Krylov subspace methods are able to solve some symmetric positive definite systems in only iterations, where is the number of unknowns in the system [6]. However, scattering problems are often described by indefinite Helmholtz equations, which are generally hard to solve using iterative methods.
This paper focuses on calculating the far field map resulting from Helmholtz and Schrödinger type scattering problems [2], which yields a representation of the scattered wave amplitude at large distances from the object of interest. The calculation of the far field map can typically be considered a twostep process. First, a Helmholtz problem with absorbing boundary conditions is solved on a finite numerical box covering the object of interest. This step generally is the main computational bottleneck in far field map computation, since discretization of a Helmholtz equation leads to a highly indefinite linear system that is notoriously hard to solve using the current generation of iterative methods. In particular, the highly efficient iterative multigrid method [8] is known to break down when applied to these type of indefinite Helmholtz problems. The observed instability of both the coarse grid correction and relaxation scheme is due to closetozero eigenvalues of the discretized operator on some intermediate multigrid levels [15]. In the second step a volume integral over the Green’s function and the numerical solution is calculated to obtain the angulardependent far field amplitude map. The twostep strategy was successfully applied to calculate impact ionization in hydrogen [35] and double photoionization in molecules [43] described by the Schrödinger equation, which in this case translates into a 6dimensional Helmholtz problem. However, computational overhead due to repeatedly solving the Helmholtz systems is significant, and supercomputing infrastructure was previously required to perform this type of calculations.
In this paper we propose a new method for the calculation of the far field map that aims at bypassing the computational bottleneck in solving the Helmholtz equation. Using basic complex analysis, the method reformulates the far field integral over the Green’s function on a complex contour. The advantage of this reformulation is that the far field integral now requires the solution of the Helmholtz equation along a complex contour, which corresponds to a damped equation, instead of requiring the realvalued scattered wave solution. The problem of solving a damped Helmholtz equation is wellknown in the literature around Helmholtz preconditioners to be significantly easier than its realvalued counterpart.
Indeed, over the past decade significant research has been performed on the construction of good preconditioners for Helmholtz problems. These results prove to be valuable in the context of this paper, albeit not in a preconditioning setting. Recent work related to this topic includes the waveray approach [9], the idea of separation of variables [33], algebraic multilevel methods [7], multigrid deflation [37] and a transformation of the Helmholtz equation into an advectiondiffusionreaction problem [20]. Moreover, in 2004 the Complex Shifted Laplacian (CSL) was proposed by Erlangga, Vuik and Oosterlee [17] as an effective preconditioner for Helmholtz problems. The key idea behind this technique is to formulate a perturbed Helmholtz problem that includes a complexvalued wavenumber. Given a sufficiently large complex shift, this implies a damping in the problem, thus making the perturbed problem solvable using multigrid in contrast to the original Helmholtz problem with realvalued wavenumbers. The concept of CSL has been further generalized in a variety of papers among which [16]. Recently a variation on the Complex Shifted Laplacian scheme by the name of Complex Stretched Grid (CSG) was proposed in [34], introducing a complexvalued grid distance instead of a complexvalued wavenumber in the perturbed system. It was furthermore shown that the CSG system is generally equivalent to the CSL scheme, and thus allows for a fast and scalable solution using classical multigrid methods.
These observations prove particularly useful in the context of far field map computation, since the damped Helmholtz equation appears naturally in the reformulation of the far field map integral proposed in this work. The level of damping is governed by the size of the complex shift (CSL) or rotation angle (CSG), which is wellknown to be bounded from below by the requirement of a stable multigrid solver, see [14]. On the other hand, choosing the shift too large may negatively impact the accuracy of the integral quadrature, hence imposing an upper bound on the damping parameter.
To validate our approach, the method is successfully illustrated on both 2D and 3D Helmholtz and Schrödinger equations for a variety of discretization levels. The absorbing boundary conditions used in this paper are based on the principle of Exterior Complex Scaling (ECS) that was introduced in the 1970’s [1], and is nowadays frequently used in scattering applications. This method is equivalent to a complex stretching implementation of Perfectly Matched Layers (PML) [5].
The outline of the article is the following. In Section 2 we introduce the notation and terminology that will be used throughout the text. Additionally, we illustrate the classical calculation of the far field map for Helmholtz type scattering problems. Section 3 contains the main theoretical insights presented in this work. Here we introduce an alternative way of calculating the far field mapping based upon a reformulation of the integral over a complex contour, for which the corresponding Helmholtz system is very efficiently solved iteratively. The new technique is validated in Section 4, and convergence results are shown for a variety of Helmholtz type model scattering problems in both two and three spatial dimensions. We demonstrate that the method allows for a very fast and scalable far field map calculation. In Section 5, the method is extensively tested on several semirealistic Schrödinger type model problems in 2D and 3D respectively. Benchmark problems include quantummechanical model problem applications for which single, double and/or triple ionization occur. Finally, along with a discussion on the topic, conclusions are drawn in Section 6.
2Helmholtz equation and far field map
In this section we introduce the general notation used throughout the text and we illustrate the classical derivation of the far field scattered wave solution and the calculation of its amplitude from a general Helmholtz type scattering problem, see [13]. The theory presented in this section provides the foundation for the work in [26], where the below results were applied to a specific quantummechanical problem (cf. Section 5).
2.1Notation and basic concepts
The Helmholtz equation is a simple mathematical representation of the physics behind a wave scattering at a certain object. Let the object be defined on a compact support area located within a domain . The equation is given by
with dimension , where is the Laplace operator, designates the right hand side or source term, and is the (spatially dependent) wavenumber, representing the material properties inside the object of interest. Indeed, the wavenumber function is defined as
where is a scalar constant denoting the wavenumber outside the object of interest. The scattered wave solution is given by the unknown function . Throughout the text we will use the following convenient notation
such that . Note that the function is trivially zero outside the object of interest where the spacedependent wavenumber is reduced to . Defining the incoming wave as , where is a unit vector that defines the direction, the righthand side is typically given by
The above expression follows directly from the fact that the total wave satisfies the homogeneous Helmholtz equation , and the incoming wave trivially satisfies . Reformulating (Equation 1), we obtain
This equation is typically formulated on a bounded domain with outgoing wave boundary conditions on , and can in principle be solved on a numerical box (i.e. a discretized subset of ) covering the support of , with absorbing boundary conditions along all edges. Let us assume that the numerical solution satisfying (Equation 2) has been calculated and is denoted by .
2.2Classical derivation of the far field map
In order to calculate the far field scattered wave pattern, the above equation is reorganized as
Note that we can replace the function in the right hand side of this equation with the numerical solution obtained from equation . In doing so, the above equation becomes an inhomogeneous Helmholtz equation with constant wavenumber
where the short notation is introduced for readability and notational convenience. It holds that for . The above equation can easily be solved analytically using the Helmholtz Green’s function , i.e.
Since the function is only nonzero inside the numerical box that was used to solve equation (Equation 2), the above integral over can be replaced by a finite integral over
In practice, this expression allows us to calculate the scattered wave solution in any point outside the numerical box, using only the information inside the numerical box.
Given the integral expression (Equation 6), the asymptotic form of the Green’s function can be used to compute the far field map of the scattered wave . In the following, this will be illustrated for a 2D model example where the Green’s function is given explicitly by
where represents the imaginary unit and is the 0th order Hankel function of the first kind. An analogous derivation can be performed in 3D, where we mention for completeness that the Green’s function is given by
To calculate the angular dependence of the far field map, the direction of the unit vector is introduced that is in 2D defined by a single angle with the positive horizontal axis, i.e. . Rewriting the spatial coordinates in polar coordinates as , the asymptotic form of the Green’s function for is given by
where we have used the fact that the Hankel function is asymptotically given by
This leads to the following asymptotic form of the 2D scattered wave solution
for . The above expression is called the 2D far field wave pattern of , with the integral being denoted as the far field (amplitude) map
The value of the integral depends only on the direction (or, in 2D, on the angle ) and the wavenumber . Expression (Equation 7) readily extends to the dimensional case, where it holds more generally that
for a function which is known explicitly and a far field map given by (Equation 8). Note that this far field map is in fact a Fourier integral of the function .
2.3Comments on the classical derivation
It is clear from the above derivation that the calculation of the far field wave pattern of a scattered wave consists of two main steps. First, one has to solve Helmholtz equation (Equation 2) with a spatially dependent wavenumber on a numerical box with absorbing boundary conditions. Once the numerical solution is obtained, it is followed by the calculation of a Fourier integral (Equation 8) over the aforementioned numerical domain. The main computational bottleneck of the calculation generally lies within the first step, since this requires an efficient and computationally inexpensive method for the solution of the indefinite Helmholtz system with absorbing boundary conditions.
The statement of the far field map presented in this section relies on the fact that the object of interest, represented by the function , is compactly supported. In particular, this is used when computing the numerical solution to equation (Equation 2) on a bounded numerical box that covers the support of . The above reasoning can however be readily extended to the more general class of analytical object functions that vanish at infinity, i.e. where . Indeed, due to the existence of smooth bump functions [23], functions with compact support can be shown to be dense within the space of functions that vanish at infinity. Consequently, every analytical function can be arbitrarily closely approximated by a series of compactly supported functions . This in turn implies that the corresponding solutions on a limited computational box can be arbitrarily close to the solution of the Helmholtz equation generated with the analytical object of interest . Intuitively, this means that if is analytical but sufficiently small everywhere outside , the computational domain may be retricted to a numerical box covering as if was compactly supported. Hence, the far field map (Equation 8) is welldefined for analytical functions that vanish at infinity. This observation will prove particularly useful in the next section.
3Far field integral: complex contour formulation
In this section we illustrate how the far field integral (Equation 8) can be reformulated as an integral over a complex contour. This new insight is consequently shown to be particularly useful with respect to the numerical computation, as it allows replacing equation (Equation 2) by a damped Helmholtz equation.
3.1Reformulation to a complex contour
The far field integral (Equation 8) can be split into a sum of two contributions: , with
Calculation of the integral is generally easy, since it only requires the expression for the incoming wave which is known analytically. The second integral however requires the solution of the Helmholtz equation on the domain (numerical box), which is notoriously hard to obtain using iterative methods. Hence, the calculation of the scattered wave forms a major bottleneck in the traditional calculation of the far field map.
However, if both and are analytical functions the integral can be calculated over a complex contour rather than the real domain as follows. Let us define a complex contour along the rotated real domain , where is a fixed rotation angle, followed by the curved segment , as presented schematically on Figure 1 for a 1D domain. The extension of the domain to multiple dimensions is straightforward, see [4]. Integral can then be written as
The second term in the above expression vanishes, since the function is per definition zero everywhere outside the object of interest , thus notably in all points . Hence, we obtain
Note that for the exponential of is increasing in all directions. At the same time the scattered wave solution , which consists of outgoing waves on the complex domain , is decaying in all directions. Additionally, the function is presumed to have a bounded support (or vanish at infinity, see Section 2.3), making the above integral computable on a limited numerical domain.
Formulation (Equation 11) of the integral is theoretically equivalent to the original formulation ( ?), since both formulations result in the same value for the integral. However, the reformulation to the complex contour provides a significant advantage from a computational point of view. Indeed, (Equation 11) indicates that the far field map can (at least partially) be computed over the full complex contour , that is, a rotation of the original realvalued domain over an angle in all spatial dimensions. Contrary to the original integral formulation ( ?), which required the scattered wave evaluated along the real domain, formulation (Equation 11) now requires the scattered wave along the complex contour. Consequently, for our new approach, the first step in the far field map calculation consists of solving the Helmholtz equation (Equation 2) on a complex contour, i.e.
which is known to be much easier to solve than the original realvalued Helmholtz equation. Classical multigrid methods allow for a fast solution of damped Helmholtz systems, see [17].
The above integral reformulation (Equation 11) implies no restrictions on the value of the rotation angle. Indeed, the rotation angle can in principle be chosen anywhere in the interval , as the equivalence between formulation ( ?) and (Equation 11) theoretically holds for any angle . However, it will be shown in Section Section 3.3 that in practice, a lower bound on the rotation angle is implied by the numerical computation of the scattered wave in the first step of the far field map computation. Note that, on the other hand, the principal incentive for not choosing excessively large is the stability of the numerical integration scheme. When the angle is very large, the scattered wave is heavily damped in all regions remote from the origin, and the integration scheme typically requires additional evaluation points near the origin, cf. [22].
3.2Solving the Helmholtz equation on a complex domain
We now show that the formulation of the Helmholtz problem on a complex rotated domain like is very similar to a complex shifted Laplacian system [17]. Indeed, both formulations can be shown to be generally equivalent, and hence are equally efficiently solvable. Consider a Helmholtz problem on a complex rotated grid of the form (Equation 12)
This equation is discretized using finite differences on a dimensional Cartesian grid with a complex valued grid distance (with ) in every spatial dimension, yielding a linear system of the form
with and where is the matrix operator expressing the stencil structure of the Laplacian. For example, discretization of the 2D Laplacian using second order finite differences yields , where the size of intrinsically depends on . Rewriting the rotation as , and multiplying both sides in by , we obtain the equivalent system
The lefthand side matrix operator of this equation is a discretization of the complex shifted Helmholtz operator . We momentarily assume to be realvalued; an analogous argument can be used for complexvalued wavenumbers under proper conditions. Hence, (Equation 15) is the discrete representation of a CSL system, which is known to be solvable using multigrid given a sufficiently large complex shift [16]. Note that in the CSL literature, the real part of the shift factor is commonly chosen as , yielding a shift of the form [17]. We expound on the relation between the complex shift and the rotation angle in Section 3.3.
The scheme given by (Equation 14) is known as Complex Stretched Grid (CSG), and was shown in [34] to hold very similar properties with respect to multigrid convergence compared to the complex shifted Laplacian system. Moreover, it is shown above that the CSL and CSG schemes are generally equivalent, and thus can be solved equally efficiently using a multigrid method.
3.3On the rotation angle
The choice of a sufficiently large complex shift parameter is vital to the stability of the multigrid solution method for problem , or equivalently . When the shift parameter is chosen below a certain minimal value denoted as , the multigrid scheme tends to be unstable and convergence will break down. A typical rule of thumb for the choice of the complex shift suggested in the CSL literature is [14]. Note that the lower limit is based on a multigrid Vcycle with standard weighted Jacobi or GaussSeidel smoothing. However, more advanced iterative techniques like ILU [40] or GMRES [11] may alternatively be used as a smoother substitute in the multigrid solver.
The requirement of a minimal shift for multigrid stability on the CSL problem can be directly translated into a minimal rotation angle for the complex scaled system. Writing the complex shift in polar notation
where and , one readily obtains
Denoting the minimal value of the complex shift by , the rotation angle must satisfy
For the rule of thumb stated above, this implies that . Note that when substituting the standard stationary multigrid relaxation scheme (Jacobi, GaussSeidel) by a more robust iterative scheme like GMRES(m) (with or ), the rotation angle is typically chosen at around , resulting in a stable multigrid scheme.
In this paper we have chosen to link the grid rotation angle to the standard ECS absorbing layer angle , as shown in Figure Figure 2. This is in no way imperious for the functionality of the method, but it appears quite naturally from the fact that both angles perturb (part of) the grid into the complex plane. Supposing the ECS boundary layer measures one quarter of the length of the entire real domain in every spatial dimension, which is a common choice, we readily derive the following relation between the rotation angle and the ECS angle :
Table 1 shows some standard values of the ECS angle and corresponding values according to (Equation 17). Following (Equation 16), we note that for a multigrid scheme with Jacobi or GaussSeidel smoothing to be stable, should be chosen no smaller than . Using the more efficient GMRES(3) method as a smoother replacement, an ECS angle around suffices to guarantee stable multigrid convergence.
On the other hand, taking too large has some drawbacks which are reflected in the accuracy of the numerical scheme for the calculation of the farfield integral. For , the solution is an exponentially decaying function towards the boundary. At the same time, the function is exponentially growing towards the boundary. Their product, which appears in the integral, remains bounded. However, when is chosen to be very large, the difference in magnitude between the two integrand factors can affect the accuracy of the numerical integral. Indeed, when multiplying two floating point numbers, one being extremely small and the other extremely large, and each is represented up to machine precision, the numerical error on this multiplication can be very large. These considerations somewhat limits the choice of the rotation angle for the complex contour.
(rad.)  

(deg.)  
4Numerical results for 2D and 3D Helmholtz problems
In this section, we validate the theoretical result presented above by a number of numerical experiments in both two and three spatial dimensions. We make use of the known efficiency of multigrid in solving damped Helmholtz equations to compute the solutions to the Helmholtz systems required in the complex valued far field integral.
The model problem used throughout this section is a Helmholtz equation of the form (Equation 2) with . The equation is discretized on a point uniform Cartesian mesh covering a square numerical domain using second order finite differences. Note that the use of a different discretization scheme would not fundamentally affect the results presented in this work. In the 2D case the spacedependent wavenumber is defined as
i.e. the object of interest takes the form of two circular pointlike objects with mass concentrated at the Cartesian coordinates and , see Figure 3 (top panel). For the 3D model problem, the following straightforward extension of the object is used
representing two spherical pointlike objects in 3D space, see Figure Figure 4 (left panel). The incoming wave scattering at the given object is defined by
where is the unit vector in the direction.

1/4  10 (10)  9 (59)  9 (560)  9 (4456)  9 (35165) 
0.24  0.20  0.21  0.20  0.20  
1/2  12 (12)  10 (63)  10 (611)  10 (4937)  9 (35405)  
0.31  0.24  0.22  0.23  0.21  
1  7 (8)  13 (83)  11 (691)  10 (4899)  10 (38975)  
0.13  0.32  0.27  0.24  0.24  
2  2 (4)  8 (54)  13 (809)  11 (5418)  10 (38051)  
0.01  0.14  0.33  0.27  0.24  
4  1 (3)  2 (17)  7 (457)  13 (6337)  11 (41848)  
0.01  0.01  0.12  0.33  0.26  
Figure 3 validates the equivalence between the classical realvalued far field map integral from Section 2 and the new complex contour formulation presented in Section 3. The 2D Helmholtz model problem with wavenumber given by (Equation 18) and is solved for using respectively a standard LU factorization method on the real domain with ECS complex boundary layers () along the domain boundary , and a series of multigrid V(1,1)cycles with Jacobi smoothing on the full complex domain () with a residual reduction tolerance of 1e6
. The standard multigrid intergrid operators used in this work are bilinear interpolation and full weighting restriction. The moduli of the wavenumber function (top) and the resulting solution (mid) are shown on Figure Figure 3 for both methods. Note how the solution on the full complex contour (right) is indeed heavily damped compared to the solution on the real domain (left). Consequently, using the numerical solution , the 2D far field map integral (Equation 8) can be calculated using any numerical integration scheme over the real or complex domain respectively. The resulting far field map is shown as a function of the angle on Figure Figure 3 (bottom). One observes that the mapping is indeed identical when calculated over the realvalued (left) and complexvalued (right) domain, conform with the theoretical results. The normalized difference between both far field computations does not exceed 0.177% (in norm), which is effectively of the same order of magnitude as the normalized error. However, the computational cost of the realdomain method for calculation of the far field map is reduced significantly by the ability to apply multigrid to the equivalent complex scaled problem.

1/4  8 (11)  6 (52)  5 (384)  5 (3190)  5 (25241) 
1.93e9  1.77e9  2.46e9  2.68e9  3.00e9  
1/2  9 (12)  8 (68)  6 (452)  6 (3392)  6 (30215)  
1.27e8  1.87e9  3.37e9  1.30e9  1.25e9  
1  5 (8)  9 (68)  8 (572)  7 (4013)  6 (30747)  
1.33e8  1.51e8  4.07e9  1.76e9  3.66e9  
2  1 (5)  5 (43)  9 (600)  8 (4456)  7 (36367)  
5.99e13  1.18e8  1.91e8  5.28e9  2.67e9  
4  1 (4)  1 (18)  5 (357)  9 (5038)  8 (39038)  
8.90e20  2.86e13  5.19e9  1.97e8  4.65e9  
In Table 2 convergence results are shown for the solution of the 3D scattered wave equation (Equation 2) using a series of multigrid V(1,1)cycles on various grid sizes. Note that the multigrid method scales perfectly as a function of the number of grid points, as doubling the number of grid points in every spatial dimension does not increase the number of Vcycles required to reach a fixed residual tolerance of 1e6
. This is a standard result from multigrid theory. Additionally and more importantly, remarkable scalability is measured for the multigrid solution method on the complex contour. Indeed, the multigrid convergence factor (and thus the corresponding work unit load required to solve the problem up to a given tolerance) is almost fully independent of the wavenumber , as can be observed from the table. From a physicalnumerical point of view it is only meaningful to consider discretizations satisfying the criterion for a minimum of 10 grid points per wavelength, cf. [4], for which the corresponding values are designated in Table 2 by a bold typesetting.
Ultimately, the computed scattered wave solution on the 3D complex domain can be used to calculate the far field integral (Equation 8). The resulting 3D far field mapping for the model problem with is shown in Figure 4. The left hand side panel shows an isosurface visualization of the 3D object of interest given by (Equation 19). On the right panel a spherical projection of the resulting 3D far field mapping is shown. The color hue indicates the value of the far field amplitude in each outgoing direction.
Note that the calculation of the scattered wave solution can be optimized even further by considering the Full Multigrid (FMG) scheme. This is a nested iteration of standard Vcycles, where on each level a series of V(1,1)cycles is used to approximately solve the error equation and supply a corrected initial guess for a finer level by interpolating the corresponding coarse grid solution.
Table 3 shows convergence results for the solution of the 3D scattered wave equation (Equation 2) using an FMG scheme. The setting is comparable to that of Table 2. A residual reduction tolerance of is imposed for each wavenumber and at every level of the FMG cycle, yielding a fine grid residual of order of magnitude . Note that the number of Vcycles performed on each level in the FMG cycle is decaying as a function of the growing grid size due to the increasingly accurate initial guess, resulting in a relatively small number of Vcycles (five to nine) to be performed on the finest level. Consequently, the number of work units (and thus the computational time) required to reach the designated residual reduction tolerance is significantly lower than the work unit load of the pure Vcycle scheme displayed in Table 2.


CPU time  0.20 s.  0.78 s.  6.24 s.  53.3 s.  462 s. 
3.3e5  7.9e5  2.7e5  1.1e5  4.6e6  
Timing and residual results from a standard FMG sweep performing only one V(1,1)cycle on each level on the 3D Helmholtz scattering problem with a moderate wavenumber are shown in Table 4 for different discretizations. Note that timings were generated using a basic nonparallelized Matlab code, using only a single thread on a simple midrange personal computer (system specifications: see caption Table 4) and taking less than 8 minutes to solve a 3D Helmholtz problem with grid points in every spatial dimension.
5Application to Schrödinger equations
This section illustrates the application of the proposed complex contour method to Schrödinger equations that are used to describe quantum mechanical scattering problems. The dimensional timeindependent Schrödinger equation for a system with unit mass is given by
where is the dimensional Laplacian, is a scalar potential, is the wave function and is the righthand side, which is often related to the ground state of the system. Depending on the total energy , the above system allows for scattering solutions, in which case the equation can be reformulated as a Helmholtz equation of the form
where the spatially dependent wavenumber is defined by . The experimental observations from this type of quantum mechanical systems are typically far field maps of the solution [43]. Indeed, in an experimental setup, detectors are commonly placed at large distances from the object compared to the size of the system. These detectors consequently measure the probability of particles escaping from the system in certain directions. In many quantum mechanical systems the potential is an analytical function, which suggests analyticity of the wavenumber in the above Helmholtz equation. Additionally, the potential is often a decaying function. Hence, for these types of problems, the wavenumber naturally satisfies all conditions for the use of the proposed complex contour method to efficiently calculate the corresponding far field map.
In the paragraphs below, we first discuss a 2D model problem where single and double ionization occur, corresponding to waves describing respectively a single particle or two particles escaping from the quantum mechanical system. The first leads to very localized evanescent waves that propagate along the boundaries of the computational domain, while the latter gives rise to waves traversing the full domain. The corresponding 2D Schrödinger problem will be solved on a discretized numerical domain for a range of energies below and above the double ionization energy threshold. For each energy, we extract the single and double ionization cross sections, which correspond to probabilities of particles escaping from the system, with the help of an integral of the Green’s function over the numerical box. The cross sections are calculated using both a traditional method, where the Helmholtz equation is solved on a standard ECSbounded grid [26], and the new complex contour method, introduced in Section 3.
The main purpose of the calculations in Sections Section 5.1 to Section 5.4 is to validate the complex contour method when applied to Schrödinger’s equation. In Section Section 5.5 the multigrid performance for solving the 2D complexvalued scattered wave system is benchmarked. Note, however, that these 2D problems essentially do not require a multigrid solver, since a direct sparse solver performs well for these relatively smallsize problems. In Section 5.6 we illustrate the convergence of a multigrid solver on a 3D Schrödinger equation, with energies that allow for triple ionization as well as double and single ionization. The previous attempts to solve these problems with the help of the complex shifted Laplacian as a preconditioner to a general Krylov method showed a notable deterioration in the convergence behavior in function of the total energy [34]. It will be shown that the new complex contour method, which allows to use multigrid as a solver, performs well for these problems.
Although the benchmark problems considered in this section mainly use model potentials, we believe that the calculations presented below are an important step towards the application of the method on realistic quantum mechanical systems.
5.1Cross sections of the 2D Schrödinger problem
Our primary aim is to validate the applicability of the new complex contour method on the 2D Schrödinger equation describing a quantum mechanical scattering problem. This problem originates from the expansion of a 6D scattering problem in spherical harmonics, see [2], in which each particle is expressed in terms of its spherical coordinates, resulting in a coupled system of 2D equations. Note that the differential operators only appear on the diagonal blocks of the system. These diagonal blocks then form the twodimensional Schrödinger equation
with boundary conditions
where represents the 2D Laplacian, and are the onebody potentials, is a twobody potential and is the total energy of the system. Since the arguments and are in fact radial coordinates in the partial wave expansion, homogeneous Dirichlet boundary conditions are implied at the and boundaries. The potentials , and are generally analytical functions that decay as the radial coordinates and become large.
Depending on the strength of the onebody potentials and , the problem allows for socalled single ionization waves, which are localized evanescent waves that propagate along the edges of the domain. We refer the reader to Sections Section 5.2 and Section 5.3 for a more detailed physical clarification. In the following, we expound on the situation with a strong attractive potential in the direction; the case with a strong potential is completely analogous. If the attraction of is strong enough, there exists a onedimensional eigenstate for every negative eigenvalue , characterized by a onedimensional Helmholtz equation
Note that and .
The far field maps of this system are then again Green’s integrals over the solution, see [26]. Indeed, the single ionization amplitude , which represents the total number of single ionized particles, is given by
where , is a onebody eigenstate, i.e. a solution of equation with a corresponding eigenvalue , and the function is a regular, normalized solution of the homogeneous Helmholtz equation
where and normalized with .
Similarly, the double ionization cross section , which measures the total number of double ionized particles, is defined by the integral
where both and are solutions to , with and , respectively, and such that . The total double ionization cross section is defined as the integral
where
The above integral expressions are obtained through a reorganization similar to the one performed in Section 2, see . For example, to calculate the single ionization cross section, equation is reorganized as
Since the left hand side is separable, the corresponding Green’s function allows us to write
Using the asymptotic form of the Green’s function, the above ultimately results in integral formulation . The double ionization integral expression can be derived in a similar way.
5.2Spectral properties of the 2D Schrödinger problem
To obtain insight in the multigrid convergence for the twodimensional Schrödinger problem, we briefly discuss the spectral properties of the discretized Schrödinger operator. The discretized 2D Hamiltonian corresponding to equation can be written as a sum of two Kronecker products and a twobody potential, i.e.
where is the onedimensional Hamiltonian, discretized using finite differences. When the twobody potential is weak relative to the onebody potentials, the eigenvalues of the 2D Hamiltonian can be approximated by
Hence, to form a better understanding of the spectral properties of , let us first consider the eigenvalues of . After discretization using second order finite differences, the onedimensional Hamiltonian can be written as a tridiagonal matrix, where the stencil
approximates the second derivatives, and the potential is a diagonal matrix evaluated in the grid points. The spectrum of closely resembles the spectrum of the Laplacian , however the presence of the potential modifies the smallest eigenvalues. The resulting spectrum is shown on the top left panel of Figure ?, which presents a closeup of the eigenvalues near the origin. A single negative eigenvalue can be observed, which is due to the attractive potential. The remaining spectrum consists of a series of positive eigenvalues located along the positive real axis. The top right panel of Figure ? shows the eigenvalues of , discretized along a complexvalued contour, i.e. the real grid rotated into the complex plane by . The grid distance used is now , which results in the following stencil for the second derivative
This implies that the spectrum of the Laplacian is rotated down into the complex plane by an angle . Figure ? shows that most of the eigenvalues of are rotated downwards over , with the exception of the bound state eigenvalue , which remains on the negative real axis. The changes of spectrum as a result of a rotation are well known in the physics literature, see for example [29].
We consequently turn to the twodimensional problem setting, where the eigenvalues of the Hamiltonian are approximately sums of the onedimensional operator eigenvalues . The resulting eigenvalues are shown on the bottom two panels of Figure ?. Again, the eigenvalues of the 2D Hamiltonian are rotated down in the complex plane when the system is discretized along a complexvalued contour. In 2D, an isolated eigenvalue appears around , and two series of eigenvalues emerge from the real axis: a first branch of eigenvalues starting at , which originates from the sum of the negative eigenvalue of the first 1D Hamiltonian combined with all the positive eigenvalues of the second 1D Hamiltonian; and a second series of eigenvalues starting at the origin, originating from the sums of the positive eigenvalues of both onedimensional Hamiltonians.
5.3Solution types: single and double ionization
The Schrödinger equation can be written shortly as , where is the Hamiltonian and the scalar is the total energy of the system. Depending on this energy , the Schrödinger system has different types of solutions. In this section we briefly expound on the physical interpretation of these solution types, using a model problem example. This section may prove less interesting to readers who are primarily interested in the computational aspects of the solution, and as such can be skipped at will.
When the energy is larger than the smallest (negative) eigenvalue of the Hamiltonian, i.e. , onebody eigenstate solutions to arise. These eigenstates can be combined into separable waves of the form
with , that form a solution to the Schrödinger problem in the region where is small and . Indeed, when is large, the potentials and are negligibly small and the resulting Schrödinger equation becomes separable in variables, where and are the solutions of the separated operators respectively. An analogous argument holds for the case when dominates and and are negligibly small, yielding evanescent waves of the form
These separable waves are solutions of the Schrödinger system for and small, and can be derived similarly to .
Single ionization  Double ionization 
Note that we can associate such a separable wave with each eigenstate of equation that corresponds to a negative eigenvalue of the Hamiltonian. These localized waves are called single ionization waves in the physics literature, since they correspond to a quantum mechanical system in which a single particle is ionized. Single ionization waves are present in the solution as soon as the energy is above the threshold. When there is a second eigenstate with negative energy, say with , an additional single ionization wave appears in the problem as soon as . We refer to the specialized literature for a detailed discussion of the ionization process, see [19].
The left panel of Figure ? shows the solution to for a total energy . The model problem under consideration fits equation , with a righthand side given by . The onebody potentials are defined by and , yielding an eigenstate of equation with energy . The twobody potential equals . The equation is discretized on a domain using 500 grid points in every spatial dimension. An ECS absorbing layer consisting of an additional 250 grid points damps the outgoing waves along the right and top edges of the domain. Note from Figure ? how the single ionization eigenstate solutions given by  appear along the edges of the domain.
When the total energy , additional asymptotic scattering solutions to the Schrödinger equation appear. If all potentials are asymptotically zero, equation boils down to a Helmholtz equation with wavenumber for and . The resulting waves are known as double ionization waves and physically correspond to the simultaneous ejection of two particles from the quantum mechanical system. In the far field, double ionization waves behave as , with an angledependent prefactor. At the same time, it is still possible to have single ionization for positive energies , since . The right panel of Figure ? shows the solution for , clearly displaying the double ionization waves in the middle of the domain and single ionization along the edges, where the two solutions coexist. Additionally, one observes that the single ionization waves oscillate faster in the  or direction than a free wave with wavenumber , since .
5.4Complex contour method for the 2D Schrödinger problem
In physical experiments, the total number of single ionized or double ionized particles is typically observed for a range of energy levels using advanced detectors. These observations are made far away from the object and effectively measure the far field amplitudes of the solutions. The outcome of this type of experiments can be predicted by calculating the single and double ionization cross sections, using the numerical solution of equation (Equation 20), see [26]. Indeed, in order to calculate these cross sections, the numerical solution to is required, which is generally hard to obtain, especially in higher spatial dimensions. However, since the potentials , and are analytical functions, the integrals for single and double ionization can be calculated along a complex contour rather than the classical real domain, in analogy to the discussion in Section 3. Consequently, one requires the scattering solution of equation on a complex contour, which is much easier to compute.
In the following, we calculate the single and double ionization cross sections for a number of energies between and using both the classical realvalued discretization and the complex contour approach proposed in this paper. The corresponding 2D scattering problems (Equation 20) are solved on a numerical domain covered by a finite difference grid consisting of 300 grid points in every spatial dimension. Additionally, an ECS absorbing boundary layer starts at and , respectively, and implements the outgoing boundary conditions, adding an additional 150 grid points in every spatial dimension. The ECS angle is . For the complex contour method, a complex scaled grid with an overall complex rotation angle is used. Solutions to for a total energy on both the classical realvalued grid and the complex contour are presented on Figure ?. Note how the solution is damped when evaluated along the complex contour.
Figure ? shows the rate of single and double ionization as a function of the total energy . The dashed and dotted lines represent the single and double ionization amplitudes calculated using the traditional realvalued method with ECS absorbing boundary conditions [26]. The solid line is the total cross section, i.e. the sum of single and double ionization, and is calculated using the optical theorem, see [31]. One observes that single ionization occurs starting from . Double ionization only occurs when , and comprises only a fraction of the single ionization cross section (cf. Figure ?). Note how the energy of the single ionized bound states rises as the total energy grows, and remains present even when . Results obtained using the complex contour approach are indicated by the and symbols on Figure ?. In this case, the Schrödinger equation is first solved on a complex contour, yielding a damped solution as shown by Figure ? (right panel), followed by the calculation of the integrals and along this complex contour. Identical results are obtained by both calculation methods, thus validating the applicability of the complex contour approach on Schrödingertype problems.
5.5Multigrid performance on the 2D Schrödinger problem
In this section, we benchmark the performance of multigrid as a solver for the 2D Schrödinger scattering problem on a complexvalued grid. It appears that the multigrid convergence rate critically depends on the value of the total energy . Indeed, Figure ? shows the convergence rate of a standard multigrid V(1,1)cycle for the 2D model problem described in the sections above. We observe that the multigrid scheme fails to converge for total energies between and . Note that this corresponds precisely to the energy range where only single ionization occurs. However, for energy levels , where both single and double ionization occur, multigrid successfully converges.
The observed convergence behavior can be explained using the spectral properties of the Hamiltonian operator introduced in Section Section 5.2. The bottom right panel of Figure ? shows the eigenvalues of , which approximate the spectrum of , discretized along a complex contour. Changing the total energy shifts the spectrum of the linear system to the left or right. For the spectrum shifts to the right, resulting in a spectrum with realvalued eigenvalues both to the left and right of the origin. This indeed implies difficulties for both the smoother and the coarse grid correction scheme. However, when , the spectrum shifts to the left, moving all eigenvalues away from zero. This results in a spectrum that is distinctly separated from the origin, corresponding to a problem more amenable to iterative solution.
From figure ?, it might appear to the reader that multigrid is not generally efficient as a Schrödinger solver due to the poor convergence in the region. However, it is important to note that multigrid is performant in the region of physical interest. Indeed, the double ionization problem requires a full twodimensional description as stated above, requiring multigrid to converge for energies . In contrast, the purely single ionized problem can be solved in a onedimensional Helmholtz setting, see , where multigrid can indeed be shown to perform well for energy levels , cf. Section 4.
Although single ionization waves are present in the solution for energy levels , they do not undermine the multigrid convergence in this regime. This is remarkable, because single ionization waves are very localized evanescent waves along the edges of the domain that generally cannot be represented efficiently on coarser grids, since there might not be enough grid points covering these regions. However, despite the fact that the coarsening strategy of the multigrid method used in this work is not adapted to evanescent waves, the damping implied by the complex contour evaluation ensures good multigrid performance.
5.6Solutions of a 3D Schrödinger equation
As demonstrated on a 2D model problem in the previous sections, the far field map (cross section) of a general Schrödinger problem can be accurately calculated using the complex contour approach. In this section we focus on the numerical solution of the 3D Schrödinger equation on the complex contour using multigrid. Note that in the threedimensional case, the use of a direct solver is strongly prohibited due to the size of the problem.
We consider the 3D Schrödinger equation, modeling a realistic scattering problem that includes single, double and triple ionization. As discussed above, this problem features very localized waves that require a sufficiently highresolution representation. The model problem derived from a 9D problem through a partial wave expansion is
with boundary conditions
We discuss a system in which , and are identical onebody potentials and , and are, similarly, identical twobody potentials. Let the strength of the onebody potential be such that there is a single negative eigenvalue for the 1D subsystem
with , where we have dropped the subscript on the 1D potential . If the twobody potential is negligibly small, then there automatically exists a bound state of the 2D subsystem. Indeed, the state is an eigenstate of the separable Hamiltonian with eigenvalue . In the presence of a small but nonnegligible twobody potential, this state will be slightly perturbed, resulting in an eigenstate that fits the 2D subsystem
The corresponding eigenvalue is . This ordering is typical for realistic atomic and molecular systems [3]. Similarly, the 3D system will have an eigenstate that looks approximately like , or any of its coordinate permutations. This 3D eigenstate fits the equation
where .
Indefinite  No  Yes  Yes  Yes  Yes 
Single ion.  No  No  Yes  Yes  Yes 
Double ion.  No  No  No  Yes  Yes 
Triple ion.  No  No  No  No  Yes 
Assuming that the potentials are such that , there are now four possible regimes of interest in equation (Equation 29), depending on the total energy . First, for , the problem is positive definite, and hence easy to solve numerically. However, in this regime no interesting physical reactions occur. Similarly, for , there are no scattering states in the solution. For energy levels , single ionization scattering occurs. Consequently, in this regime, there exist scattering solutions that are localized along one of the three axes in the 3D domain. These solutions take the form as , where is the eigenstate of and is a scattering solution satisfying outgoing wave boundary conditions. Similar solutions are found for the respective coordinate permutations. For energies , both single and double ionization occurs. The solution contains — besides single ionization waves — double ionization waves of the form , where is an eigenstate of and is a 2D scattering state satisfying the outgoing wave boundary conditions. Together with coordinate permutations, these waves are localized along the faces of the 3D domain, where one of the three coordinates, , or , is small. Finally, for , the solution contains, in addition, triple ionization waves. These are waves that describe a quantum mechanical system that is fully broken up into its subparticles. In this case, all three relative coordinates , and can become large, resulting in a wave that extends to the entire domain.
Note that in fact only the latter problem, when and triple ionization is present, requires a full 3D description. For the regimes in which only single ionization occurs, a description in the form of a coupled set of 1D equations is sufficient, due to the separated character of the solution. Similarly, for problems with both single and double ionization, but no triple ionization, a simpler 2D description such as the one given by can be used to fully describe the physics behind the problem. Hence, in view of the efficient solution of the 3D Schrödinger problem , our main interest goes out to the regime.
5.7Multigrid performance on the 3D Schrödinger problem
We now study the convergence of a multigrid solver for the 3D Schrödinger equation (Equation 29) for energies that cover all possible scattering regimes. The model problem under consideration is a straightforward generalization of the 2D model presented in Section 5.3, featuring onebody potentials , and