Application of the spectral element method to the solution of the multichannel Schr\"odinger equation

TL;DR

This paper demonstrates the application of the spectral element method to accurately solve the multichannel Schrödinger equation, efficiently computing bound states and scattering matrices with fewer grid points.

Contribution

The paper introduces a spectral element approach for multichannel Schrödinger equations, offering high accuracy and computational efficiency over traditional methods.

Findings

Spectral accuracy achieved with fewer grid points.

Efficient solution of scattering matrices using sparse matrix solvers.

Comparison shows advantages over log-derivative propagators.

Abstract

We apply the spectral element method to the determination of scattering and bound states of the multichannel Schr\"odinger equation. In our approach the reaction coordinate is discretized on a grid of points whereas the internal coordinates are described by either purely diabatic or locally diabatic (diabatic-by-sector) bases. Bound levels and scattering matrix elements are determined with spectral accuracy using relatively small numbers of points. The scattering problem is cast as a linear system solved using state-of-the-art sparse matrix non iterative packages. Boundary conditions can be imposed so to compute a single column of the matrix solution. A comparison with log-derivative propagators customarily used in molecular physics is performed. The same discretization scheme can also be applied to bound levels that are computed using direct scalable sparse-matrix solvers.

Tables

Table 1. Table 1: CPU (in seconds) for various P 𝑃 P and h ℎ h combinations corresponding to two given values of the error on δ 𝛿 \delta extracted from fig.  3 . The columns E c ​ o ​ l subscript 𝐸 𝑐 𝑜 𝑙 E_{col} ind. and E c ​ o ​ l subscript 𝐸 𝑐 𝑜 𝑙 E_{col} dep. correspond to the CPU time of steps to be done once per collisional energy ( E c ​ o ​ l subscript 𝐸 𝑐 𝑜 𝑙 E_{col} ind.) and at each collision energy ( E c ​ o ​ l subscript 𝐸 𝑐 𝑜 𝑙 E_{col} dep.).

	${\log }_{10}\left(\Delta \delta \right)=-2$			${\log }_{10}\left(\Delta \delta \right)=-6$
type	$h$	$E_{col}$ ind.	$E_{col}$ dep.	$h$	$E_{col}$ ind.	$E_{col}$ dep.
$P=3$	0.10	368	2267
$P=4$	0.38	140	795
$P=5$	0.77	95	492	0.24	282	1556
$P=6$	1.23	73	373	0.49	180	814
$P=7$	1.81	59	339	0.80	141	778
Johnson	0.10	354	80	0.01	3422	786

Equations

$$\left[-\frac{{\hslash }^2}{2\mu }\frac{{\partial }^2}{\partial R^2}+V(R,\Omega )\right]\Psi (R,\Omega )=E\Psi (R,\Omega ),$$

$${\int }_{R_1^{(m)}}^{R_{P_m}^{(m)}}f(R)dR=\sum _{p=1}^{P_m}{w}_p^{(m)}f(R_p^{(m)}),$$

$$C_i^{(m)}(R)=c_i\left(2\frac{R-R_1^{(m)}}{R_{P_m}^{(m)}-R_1^{(m)}}-1\right).$$

$$c_i(x)=\frac{-(1-x^2)}{N(N-1)P_{N-1}(x_i)(x-x_i)}\frac{d{P}_{N-1}(x)}{dx},$$

$${\mathcal{C}}_i^{(m)}(R)=C_i^{(m)}(R).$$

$${\cal C}_{1}^{(m)}(R)=\left\{\begin{array}[]{ll}C_{P_{m-1}}^{(m-1)}(R)&R\in[R_{1}^{(m-1)},R_{P_{m-1}}^{(m-1)}]\\ C_{1}^{(m)}(R)&R\in[R_{1}^{(m)},R_{P_{m}}^{(m)}]\end{array}\right.$$

$$\displaystyle\omega_{a}=\left\{\begin{array}[]{ll}\left(w_{P_{m-1}}^{(m-1)}+w_{1}^{(m)}\right)&\text{if $R_{a}$ is inter-element}\\ w_{p}^{(m)}&\text{otherwise}.\end{array}\right.$$

$$\displaystyle{\cal C}_{a}(R)=\left\{\begin{array}[]{ll}{\cal C}_{1}^{(m)}&\text{if $R_{a}$ is inter-element}\\ {\cal C}_{p}^{(m)}&\text{otherwise}.\end{array}\right.$$

$${\int }_{R_1}^{R_L}{\mathcal{C}}_b(R){\mathcal{C}}_a(R)dR={\delta }_{ba}{\omega }_a$$

$$\Psi (R,\Omega )=\sum _{a=1}^L{\Phi }_a(\Omega ){\mathcal{C}}_a(R),$$

$${\int }_{R_1}^{R_L}{\mathcal{C}}_a(R)\frac{{\partial }^2\Psi (R,\Omega )}{\partial R^2}dR$$

$${\mathcal{C}}_a(R_{P_m}^{(m)})\frac{\partial \Psi (R_{P_m}^{(m)},\Omega )}{\partial R}={\mathcal{C}}_a(R_1^{(m+1)})\frac{\partial \Psi (R_1^{(m+1)},\Omega )}{\partial R}$$

$${\int }_{R_1}^{R_L}{\mathcal{C}}_a(R)\frac{{\partial }^2\Psi (R,\Omega )}{\partial R^2}dR$$

$${\mathcal{T}}_{ab}=\sum _{m=1}^M{\int }_{R_1^{(m)}}^{R_{P_m}^{(m)}}\frac{d{\mathcal{C}}_a(R)}{dR}\frac{d{\mathcal{C}}_b(R)}{dR}dR$$

$${\mathcal{T}}_{ab}={\int }_{R_1}^{R_L}\frac{d{\mathcal{C}}_a(R)}{dR}\frac{d{\mathcal{C}}_b(R)}{dR}dR.$$

$$\sum _{b=1}^L{\mathcal{T}}_{ab}{\Phi }_b(\Omega )+\frac{2\mu }{{\hslash }^2}{\omega }_a\left[V(R_a,\Omega )-E\right]{\Phi }_a(\Omega )={\delta }_{La}{\partial }_R\Psi {|}_{R=R_L}-{\delta }_{1a}{\partial }_R\Psi {|}_{R=R_1}.$$

$${\Phi }_a(\Omega )=\sum _{\alpha =1}^{N^{(a)}}{F}_{a\alpha }{\varphi }_{\alpha }^{(a)}(\Omega ),$$

$$\sum _{b=1}^L\sum _{\beta =1}^{N^{(b)}}{\mathcal{T}}_{ab}{\mathcal{O}}_{a\alpha ,b\beta }{F}_{b\beta }+\frac{2\mu }{{\hslash }^2}\sum _{\beta =1}^{N^{(a)}}{\omega }_a\left[{\mathcal{U}}_{\alpha \beta }(R_a)-E{\delta }_{\alpha \beta }\right]F_{a\beta }={\delta }_{La}{\xi }_{\alpha L}-{\delta }_{1a}{\xi }_{\alpha 1}$$

$${\xi }_{a\alpha }=⟨{\varphi }_{\alpha }^{(a)}|{\partial }_R\Psi {|}_{R=R_a}{⟩}_{\Omega }.$$

$$\sum _{b=2}^{L-1}\sum _{\beta =1}^{N^{(b)}}{\mathcal{T}}_{ab}{\mathcal{O}}_{a\alpha ,b\beta }{F}_{b\beta }^{(n)}+\frac{2\mu }{{\hslash }^2}\sum _{\beta =1}^{N^{(a)}}{\omega }_a{\mathcal{U}}_{\alpha \beta }(R_a)F_{a\beta }^{(n)}=E_n{\omega }_a{F}_{a\alpha }^{(n)},a=2,\dots ,L-1.$$

$$\sum _{b=2}^L\sum _{\beta =1}^{N^{(b)}}{\mathcal{T}}_{ab}{\mathcal{O}}_{a\alpha ,b\beta }{\bar{F}}_{b\beta }^{(\gamma )}+\frac{2\mu }{{\hslash }^2}\sum _{\beta =1}^{N^{(a)}}{\omega }_a\left[{\mathcal{U}}_{\alpha \beta }(R_a)-E{\delta }_{\alpha \beta }\right]{\bar{F}}_{a\beta }^{(\gamma )}={\delta }_{La}{\delta }_{\alpha \gamma },a=2,\dots ,L.$$

$${\Psi }^{(I)}(R,\Omega )=\sum _{\alpha =1}^{N^{(L)}}\left[{\tilde{f}}_{\alpha }(R){\delta }_{\alpha I}-{\tilde{g}}_{\alpha }(R)K_{\alpha I}\right]{\Phi }_{\alpha }^{(L+1)}(\Omega ).$$

$${\Psi }^{(I)}(R_L,\Omega )=\sum _{\gamma =1}^{N^{(L)}}{\bar{\Psi }}^{(\gamma )}(R_L,\Omega )N_{\gamma I}$$

$${\partial }_R{\Psi }^{(I)}(R_L,\Omega )=\sum _{\gamma =1}^{N^{(L)}}{\partial }_R{\bar{\Psi }}^{(\gamma )}(R_L,\Omega )N_{\gamma I}.$$

$$\mathbf{K}={\left(\mathbf{g}-\mathbf{R}{\mathbf{g}}^{\prime }\right)}^{-1}\left(\mathbf{f}-\mathbf{R}{\mathbf{f}}^{\prime }\right).$$

$$f_{\alpha \beta }={\tilde{f}}_{\alpha }(R_L){\mathcal{O}}_{L\alpha ,(L+1)\beta },f_{\alpha \beta }^{\prime }={\tilde{f}}_{\alpha }^{\prime }(R_L){\mathcal{O}}_{L\alpha ,(L+1)\beta }$$

$$F_{L\alpha }^{(I)}=\sum _{\beta }\left[f_{\alpha \beta }{\delta }_{\beta I}-g_{\alpha \beta }{K}_{\beta I}\right].$$

$${\xi }_{L\alpha }^{(I)}=\sum _{\beta }\left[f_{\alpha \beta }^{\prime }{\delta }_{\beta I}-g_{\alpha \beta }^{\prime }{K}_{\beta I}\right].$$

$$M_{a\alpha ,b\beta }={\mathcal{T}}_{ab}{\mathcal{O}}_{a\alpha ,b\beta }+\frac{2\mu }{{\hslash }^2}{\omega }_a\left[{\mathcal{U}}_{\alpha \beta }(R_a)-E{\delta }_{\alpha \beta }\right]{\delta }_{ab}$$

$$\sum _{b=2}^{L-1}\sum _{\beta =1}^{N^{(b)}}{M}_{a\alpha ,b\beta }{F}_{b\beta }^{(I)}+\sum _{\beta =1}^{N^{(L)}}\left[M_{a\alpha ,L\beta }+{\delta }_{La}{Y}_{\alpha \beta }^g\right]X_{\beta }^{(I)}=-\sum _{\beta =1}^{N^{(L)}}{M}_{a\alpha ,L\beta }{f}_{\beta I}+{\delta }_{La}{f}_{\alpha I}^{\prime }.$$

Full text

I Introduction

The solution of partial differential equations, ubiquitous in all areas of physics, can be tackled by a variety of numerical methods developed over last decades. Solution algorithms can essentially be divided into propagation and basis expansion approaches. The former build the solution iteratively from a known initial value up to the final propagation distance, where suitable boundary conditions are imposed. Such methods are easy to implement and cheap in memory storage but provide relatively low convergence rates as a function of the step size. Large number of steps may therefore be needed, such that accumulation of roundoff error can limit the accuracy in particular for complex systems. Due to its pivotal role in quantum dynamics, the time-independent Schroödinger equations has been granted a particular attention in the molecular physics community. Popular propagation algorithms include the log-derivative propagator of Johnson Johnson (1973), later improved by Manolopoulos Manolopoulos (1986), and the renormalized Numerov algorithm Johnson (1978).

In basis expansions the solution is determined altogether as a development over a basis, usually trigonometric or polynomial, with suitable conditions imposed at the boundaries. A main advantage is the exponential numerical convergence as a function of the expansion order Boyd (2001). Grid-discretization methods are particular basis expansions in which each basis functions is nonzero at a unique grid point. In global approaches the whole interval of interest is represented as a discrete grid of points. Global grid techniques have been introduced in molecular physics in the context of the so-called discrete variable representation (DVR) Lill et al. (1982). One drawback of global grids is the need to introduce nonlinear coordinate transformations to efficiently represent complex solutions varying on disparate length scales such as ultracold processes Tiesinga et al. (1998); Kokoouline et al. (1999).

Local approaches subdivide the interval of interest in subintervals, often termed elements. Basis functions localized in subintervals are used to expand the solution. As compared to global approaches, since the resulting matrix is sparse one can apply performant sparse linear algebra packages to carry out the operations needed in the specific problem. Moreover, at least in one dimensional problems the element size can be tailored to the solution in a straightforward way.

One widespread local approach is the finite element method Hughes (2000). Being based on low order polynomial expansions, such method requires large number of points to achieve high accuracy. Use of high order polynomials as basis functions in each element marked the birth of spectral element approaches few decades ago Orszag (1980). The spectral element method is nowadays a well established tool in scientific and engineering computations Karniadakis and Sherwin (2005); Canuto et al. (2007).

Use of this computational technique in molecular physics has been pioneered for one dimensional problems in Ref. Manolopoulos and Wyatt (1988). More recently, a spectral element approach in two spatial dimensions has been presented in Rescigno and McCurdy (2000); Tao et al. (2009), where it was termed finite-element DVR. Appropriate scattering boundary conditions were imposed using a spatial rotation in the complex plane known as exterior complex scaling approach Rescigno and McCurdy (2000); Scrinzi (2010). A combination of the spectral element and the slow variable representation Tolstikhin et al. (1996) has been proposed in order to compute weakly bound states of triatomic systems in the hyperspherical framework Suno (2011). Finally, the finite-element DVR has been used as a time-independent representation in multidimensional time-dependent calculations; See e.g. Schneider et al. (2006); Telnov and Chu (2007); Hu (2010).

Main aim of the present work is to explore the usefulness of the spectral element method in quantum dynamics for time-independent multichannel problems. More particularly, we point out that combining the spectral element method with traditional molecular basis or diabatic-by-sector expansions Launay and Le Dourneuf (1989) to treat the internal coordinates optimizes sparsity and size of the discretized Hamiltonian. One major advantage is that the wavefunction is obtained at no extra computational cost. Moreover, the spectral nature of the method allows the accuracy of the solution to be estimated reliably in each region of space. Subsequent step refinements lead to a grid tailored to the interaction strength in various regions of space. We show that the spectral element formulation lends itself to imposing in a natural way different boundary conditions for scattering and bound-state calculations.

The paper is organized as follows: Sec. II presents the discretization scheme and introduces various boundary conditions, Sec. III discusses numerical results on a realistic ro-vibrational system, Sec. IV summarizes and concludes this work.

II Hamiltonian discretization

Fundamentals of the spectral-element approach can be found in textbooks and articles Canuto et al. (2007); Manolopoulos and Wyatt (1988); Rescigno and McCurdy (2000). In order to set the notation and to illustrate the specific approach we follow to combine full or locally diabatic expansions and the grid basis, we reproduce in this section the main steps of the derivation from scratch.

We consider a generic time-independent multidimensional problem and identify an “external”reaction coordinate $R$ describing the “size” of the system and a set of “internal”variables denoted collectively as $\Omega$. For instance, in the two-body problem $R$ typically represents the distance between the particles and $\Omega$ the orientation of the inter-particle position vector. In three-body systems $R$ may represent the hyperradius and $\Omega$ a set of hyperangles. Note that in general $\Omega$ may comprise coordinates with physical dimension of length, like in the case of our test atom-molecule ro-vibrational problem described in detail in Sec. III.

The time-independent Schrödinger equation is schematically written

[TABLE]

to be solved in the hyper-region $R_{\text{min}}\leq R\leq R_{\text{max}}$. Here $V$ contains various potential energy terms and/or differential operators acting on the internal variables $\Omega$. The derivation turns out to be formally simpler if, in the spirit of the slow variable representation Tolstikhin et al. (1996), $R$ is discretized first and the internal variables $\Omega$ at a second stage.

The radial interval is partitioned into $M$ subintervals or elements. We generate in each subinterval $m$ a set of $P_{m}$ Gauss-Lobatto points and weights $(R_{p}^{(m)},w_{p}^{(m)}),p=1,\dots,P_{m}$, with $R_{1}^{(m)}$ and $R_{P_{m}}^{(m)}$ the subinterval endpoints Boyd (2001). Note that since the endpoints of contiguous intervals are such that $R_{P_{m-1}}^{(m-1)}=R_{1}^{(m)}$, the number of distinct points in the complete grid is $L=\sum_{m=1}^{M}(P_{m}-1)+1$. The local points and weights can be used to implement the Gauss-Lobatto quadrature

[TABLE]

an integration rule exact for polynomials up to degree $2P_{m}-3$. Each point can be associated with a Gauss-Lobatto cardinal or shape function defined such that $C_{i}^{(m)}(R_{p}^{(m)})=\delta_{ip}$ at the nodal points inside the element and continued as identically zero outside the element, $C_{i}^{(m)}(R)=0$ if $R\notin[R_{1}^{(m)},R_{P_{m}}^{(m)}]$. The $C_{i}^{(m)}$ functions can be obtained by linear mapping in terms of the corresponding cardinal functions $c_{i}(x)$ of the primitive interval $x\in[-1,1]$

[TABLE]

Explicitly, the $c_{i}$ for $N$ grid points can be expressed in terms of the derivatives of the Legendre polynomial of order $N-1$ as follows Boyd (2001)

[TABLE]

where the Gauss-Lobatto points $x_{i}$ in the primitive interval comprise the endpoints $\pm 1$ and the $(N-2)$ zeros of the $\frac{dP_{N-1}(x)}{dx}$ polynomial.

The Gauss-Lobatto cardinal functions associated to the internal points $p=2,\cdots,P_{m}-1$ vanish at the element endpoints $R=R_{1}^{(m)},R_{P_{m}}^{(m)}$ and following Ref. Canuto et al. (2007) will be referred to as “internal functions”. We also conventionally consider as internal the cardinal functions relative to the first $R_{1}^{(1)}$ and last $R_{P_{M}}^{(M)}$ grid points. For each internal point, that is for $m=2,\ldots,M-1$ and $i=2,\ldots,P_{m}-1$ as well as for $(i,m)=(1,1)$ and $(i,m)=(P_{M},M)$, we will simply take as basis functions the cardinal functions

[TABLE]

The construction of the basis functions associated to the remaining $(M-1)$ inter-element points $R_{P_{m-1}}^{(m-1)}=R_{1}^{(m)}$ with $m=2,\cdots,M$ are obtained by “glueing”cardinal functions Canuto et al. (2007). These interface or bridge functions are defined by

[TABLE]

for $m=2,\ldots,M$. Note that contiguous subintervals are only connected through such interface functions.

In order to build a global representation of the Hamiltonian, we now introduce a single index $a=1,\ldots,L$ running over the $L$ distinct points of the full grid and note $R_{a}$ such distinct grid points. We define global weights

[TABLE]

Similarly, we build a global grid basis comprising internal and interface functions

[TABLE]

Equations. (5) and (12) guarantee that the orthogonality relation

[TABLE]

holds with at least Gaussian quadrature accuracy; See Eq. (2).

For each value of the internal coordinates $\Omega$ we now develop the system wavefunction on the radial basis

[TABLE]

the coefficients being equal to the wavefunction evaluated at the grid points $\Phi_{a}(\Omega)=\Psi(R_{a},\Omega)$. The Schrödinger equation is now projected on the basis functions ${\cal C}_{a}$. The second derivative arising from the radial kinetic energy term gives rise to an integral in $R$ that is further developed as a sum of integrals restricted to each element

[TABLE]

where one integration by parts has been performed for the second equality. Noticing that for two consecutive elements one has

[TABLE]

Eq. (II) reduces to

[TABLE]

Note that the boundary terms cancellation of Eq. (16) holds for the exact solution but is only approximatly valid when the solution is computed as an expansion on a finite basis. In other terms, the numerical solution will in general have a discontinous derivative at the element interfaces. However, such left-right discontinuity tends to zero exponentially for a converged calculation and as such does not affect the fast convergence rate demonstrated in Sec. III. Using the decomposition Eq. (14) for the evaluation of $\partial\Psi(R,\Omega)/\partial R$ one gets a term involving the matrix

[TABLE]

formally recast as

[TABLE]

Using the definitions Eq. (5) and Eq. (6), the linear mapping in Eq. (3), and approximating the integrals on the rhs by the quadrature of Eq. (2), the kinetic matrix can be expressed in terms of analytically known Gauss-Lobatto derivation matrices $C_{j}^{\prime}(x_{i})$; See e.g. Boyd (2001).

A closer look at Eq. (18) keeping into account the local character of the ${\cal C}_{a}(R)$ functions shows that most elements of ${\cal T}$ are zero. More specifically, ${\cal T}_{ab}=0$ if ${\cal C}_{a}$ and ${\cal C}_{b}$ are both internal functions and do not belong to the same element. If ${\cal C}_{a}$ is an interface function, thus at the interface of two elements $m$ and $m+1$, ${\cal T}_{ab}=0$ if $b$ does not belong to any of the two $m$ and $m+1$ elements, while ${\cal T}_{ab}\neq 0$ if ${\cal C}_{b}$ is an internal or an interface function of element $m$ or $m+1$. In addition, ${\cal T}_{ab}\neq 0$ if both $a$ and $b$ belong to the same element. The potential energy is approximately diagonal in the grid basis with diagonal elements $\omega_{a}{V}(R_{a},\Omega)$ if quadrature Eq. (2) is used. Collecting all the terms, the matrix form of the Schrödinger equation finally reads

[TABLE]

We now introduce an internal coordinate basis $\phi_{\alpha}(\Omega)$ whose nature or dimension does not need for the moment to be specified. We write therefore

[TABLE]

where the superscript $(a)$ stresses the possible dependence of the basis on the grid point. Insertion of Eq. (21) in Eq. (20) leads to the algebraic equations

[TABLE]

where ${\cal O}_{a\alpha,b\beta}=\langle\phi_{\alpha}^{(a)}|\phi_{\beta}^{(b)}\rangle_{\Omega}$ is the overlap matrix element over the $\Omega$ coordinates and ${\cal U}_{\alpha\beta}(R_{a})=\langle\Phi_{\alpha}^{(a)}(\Omega)|V(R_{a},\Omega)|\Phi_{\beta}^{(a)}(\Omega)\rangle$. The quantities $\xi_{\alpha a}$ with $a=1,L$ are the normal derivatives of the wavefunction at the integration boundaries in channel $\alpha$, i.e.

[TABLE]

Equations (22), supplemented by the proper boundary conditions in Sec. II.1 and II.2, represent the key formal result of the paper. In order to maximize sparsity one can require that the basis $\phi_{\alpha}^{(a)}$ does not depend on the grid point $a$. This is for example the case when using spherical or hyperspherical harmonics, or ro-vibrational molecular states. In this case, the ${\cal O}$ matrix reduces to the identity matrix.

A pictorial representation of the resulting Hamiltonian matrix is given in Fig. 1 where each small block correspond to fixed grid indices and varying channel indices. Similar matrix structure representations can be found elsewhere in the literature, for instance in Scrinzi (2010); Hu (2010). For this example, three subintervals $M=3$ are considered, the number of Gauss Lobatto points and the number of basis for the internal coordinates are identical for the three subelements and are fixed to $P_{m}=4$ and $N^{(a)}=5$ channels. The non-zero elements arising from the kinetic part are depicted in gray while the ones resulting from the potential are in blue. For the illustrative case presented, 340 over the 50x50 matrix elements are non-zero, which amounts to a filling factor of 13.6%.

Restoring some flexibility in the choice of the basis used for the $\Omega$ part but still leaving quite a large sparsity in the full Hamiltonian matrix can also be obtained by imposing that the basis does not vary within each element. This means that the dependence of $\phi_{\alpha}^{(a)}$ on $(a)$ is replaced by an ensemble of basis functions in which the same basis is used for all points within the same element $m$ with the exception of the inter-element points for which alternative basis functions may be used. One such example would be the set of eigenvectors obtained through diagonalization of the reduced Hamiltonian ${V}(R_{x},\Omega)$ at a fixed point $R_{x}$ inside the element (diabatic-by-sector method Launay and Le Dourneuf (1989)). The diabatic-by-sector approach trades loss of some sparsity with a (possible) reduction in the basis size.

As an example, the matrix structure for a more flexible basis choice is depicted in Fig. 2 in which different numbers of Lobatto points are used $P_{1}=4$, $P_{2}=5$ and $P_{3}=7$. In addition the first element contains 5 basis functions whereas the other two 4 alternative basis functions. The basis associated to the first inter-element point is taken to be identical to the 5 basis functions of the first element. The modification of the basis between elements 1 and 2 induces additional non-zero elements (pink in the figure) due to basis overlaps.

In the extreme case where the Hamiltonian is diagonalized at each point, one retrieves the slow-variable formulation proposed for hyperspherical bound states in Ref. Suno (2011). Such fully adiabatic procedure does optimize the basis size, but results in full overlap matrices at all off-diagonal grid elements, putting severe memory constraints on the size of treatable systems.

Realistic calculations usually require large number of elements and points to be converged. Therefore, in the locally diabatic formulation the empty part of the matrix becomes large and the filling factor decreases significantly. More quantitatively, let us consider a potential represented by a full matrix with $N$ channels. For a partition composed of $M$ elements with $P$ points per element the number of nonzero elements is $N(N+1)\left[M(P-1)+1\right]/2+P(P-1)MN/2$, that reduces to $\approx LN\left[N+P+1\right]/2$ for $M,P\gg 1$. Note that due to matrix symmetry only elements above (or below) the main diagonal have been taken into account. The total number of elements (now both above and below the main diagonal) is $N^{2}\left[M(P-1)+1\right]^{2}\approx(NLP)^{2}$, resulting in a filling factor $\approx(N+P+1)/(2NLP^{2})$ that may easily drop below 1%. We remark that this worst case scenario of a full potential matrix seldom happens in molecular physics due to the tensor nature of at least part of the interaction and to the accompanying selection rules.

Depending on the problem, purely diabatic and diabatic-by-sector representation can also be conveniently combined. For instance, in Ref. Simoni et al. (2015) the present algorithm was used to join a purely diabatic spherical harmonics basis at short range with a diabatic-by-sector representation at long range to treat the difficult problem of ultracold collisions between two polar molecules in an optical tube.

II.1 Bound states

Bound states $\Psi^{(n)}$ of the system with energy $E_{n}$ are calculated imposing that the solution of the Schrödinger equation vanishes at the boundaries $R_{1}$ and $R_{L}$ of the radial interval. Accordingly, for $a=1$ and $L$ the expansion coefficients $F_{a\alpha}^{(n)}=0$, and all elements of the discretized Hamiltonian with grid indices $a,b=1,L$ in the system of equations (22) can be dropped.

The resulting equations for the multichannel bound-state solution at the remaining $L-2$ points present themselves in the form of a generalized eigenvalue problem

[TABLE]

The $\omega_{a}$ factor on the rhs can be removed by redefining as new unknown $\sqrt{\omega_{a}}F_{a\alpha}^{(n)}$ and by right multiplication by the diagonal matrix with element $1/\sqrt{\omega_{a}}$. In this way, the problem is expressed as an ordinary eigenvalue problem. Finally note that, if needed, equations (22) with left grid index $a=1$ and $L$ can be used to compute the normal derivative components $\xi_{a\alpha}^{(n)}$ as a function of the $F_{b\alpha}^{(n)}$ with $b=2,\dots,L-1$.

II.2 Scattering states

In the case of scattering solutions, a number $N^{(L)}$ of linearly independent solutions with energy $E$ can by built from Eq. (22). We consider for definiteness the most common case where at the left end point $R_{1}$ the wavefunction vanishes. As in the bound-state problem, this implies that all lines and columns in the system of equations (22) with grid index $a=1$ can be dropped. At the other edge of the grid, we discuss below two approaches to obtain the physical wavefunction and thus the relevant scattering observables. The first one requires to compute a matrix comprising all linearly independent solutions of the discretized Schrödinger equation, the second one can be used to determine a single column-vector solution.

II.2.1 $R$-matrix boundary conditions

The so-called $R$-matrix solutions ${\bar{\Psi}}^{(\gamma)}$ are defined by the condition that their normal derivative vanishes on the surface $R=R_{L}$ in all but channel $\gamma$, where it is unity. Therefore, such $N^{(L)}$ independent $R$-matrix solutions with energy $E$ can be determined by imposing $\xi_{\alpha L}=\delta_{\alpha\gamma}$ and solving the linear system

[TABLE]

The $R$-matrix $\mathbf{R}$ is simply defined as the matrix with elements $R_{\alpha\beta}={\bar{F}}_{L\alpha}^{(\beta)}$.

Solutions with physical boundary conditions can be written for $R\geq R_{L}$ as a superposition of solutions of the asymptotic Hamiltonian, comprising angular channel eigenfunctions $\Phi_{\alpha}^{(L+1)}$ and of reference regular $\tilde{f}$ and irregular $\tilde{g}$ radial functions

[TABLE]

The channel eigenfunctions are $R$-independent and the superscript $(L+1)$ is merely introduced as an additional artificial grid point for ease of notation in subsequent formal manipulations. A solution $\Psi^{(I)}$ corresponds to a wave incoming in channel $I$ with scattered waves in all channels $\alpha$, with amplitudes $K_{\alpha I}$. If $f$ and $g$ are real standing waves the coefficients $K_{\alpha I}$ form the reaction matrix $\mathbf{K}$.

The solutions $\Psi^{(I)}$ and its normal derivative can be expressed on the surface $R=R_{L}$ as linear combinations of the $R$-matrix solutions ${\bar{\Psi}}^{(\gamma)}$ with constant coefficients $N_{\gamma I}$

[TABLE]

and

[TABLE]

Following the standard asymptotic matching procedure Launay and Le Dourneuf (1989), the expression Eq. (26) is inserted on the lhs of Eqs. (27) and (28) and the latter are projected on the angular basis $\Phi^{N^{(L)}}(\Omega)$. The resulting linear system can be easily solved for $\mathbf{K}$ in terms of $\mathbf{R}$

[TABLE]

Here, matrices $\mathbf{f}$ and ${\mathbf{f}}^{\prime}$ are respectively defined as

[TABLE]

as a function of the overlap between the asymptotic channels and the angular basis at last grid point. A similar definition holds for $\mathbf{g}$ and ${\mathbf{g}}^{\prime}$.

II.2.2 Scattering boundary conditions

Rather than going through the determination of $N^{(L)}$ independent $R$-matrix solutions, scattering boundary conditions can also be incorporated directly in the linear system of equation (22). To this aim, we first impose that at last grid point $a=L$ the wavefunction takes the form (26)

[TABLE]

Similarly, the normal derivative channel components on the rhs of Eq. (22) becomes

[TABLE]

For notational ease, we define the matrix on the lhs of Eq. (22)

[TABLE]

As it will be clear from equation (34) below, in order to obtain a symmetric linear system, it is necessary to introduce the new unknown ${\mathbf{X}}={\mathbf{g}}{\mathbf{K}}$ in the place of $\mathbf{K}$. We also define the log-derivative ratio ${\mathbf{Y}}^{g}={\mathbf{g}}^{\prime}{\mathbf{g}}^{-1}$, such that the quantity ${\mathbf{g}}^{\prime}{\mathbf{K}}$ on the rhs of Eq. (32) becomes ${\mathbf{g}}^{\prime}{\mathbf{K}}={\mathbf{Y}}^{g}{\mathbf{X}}$.

With these definitions, simple matrix algebra allows one to cast the system (22) into the form

[TABLE]

As a final step, the $K$-matrix can be computed from the definition of ${\mathbf{X}}$ by solving the linear system ${\mathbf{g}}{\mathbf{K}}={\mathbf{X}}$. It is important to notice that for a given incoming wave labeled by index $I$ one can determine a single column of the matrix solution ${\mathbf{X}}$ and thus of $\mathbf{K}$. If one uses complex algebra and replaces $f_{\alpha}$ and $g_{\alpha}$ by travelling waves $h_{\alpha}^{(-)}$ and $h_{\alpha}^{(+)}$, the asymptotic condition (26) becomes

[TABLE]

with $\mathbf{S}$ the scattering matrix, whose elements are directly related to observables. The procedure to determine $\mathbf{K}$ presented in this section applies as is to the determination of $\mathbf{S}$, leading to the equivalent of Eq. (34) with ${\mathbf{f}}$ and ${\mathbf{g}}$ replaced by ${\mathbf{h}^{(-)}}$ and ${\mathbf{h}^{(+)}}$ and ${\mathbf{Y}}^{g}$ by ${\mathbf{Y}}^{h^{(+)}}={\mathbf{h}}^{(+)\prime}[{\mathbf{h}}^{(+)}]^{-1}$. Determining a single column of interest of the scattering matrix may lead to computational advantages, in particular in problems with large numbers of open channels.

II.3 Spectral log-derivative propagation

In spite of the sparse character of the discretized Hamiltonian, memory can become a limiting factor for systems described by large numbers of collision channels. In this case, it may be necessary to split the full propagation interval in smaller intervals, each comprising for instance only one element. The scattering equation is solved in any given element element to determine at each point a matrix of linearly independent solutions ${\mathbf{F}}_{a}$ with elements $F^{(I)}_{\alpha a}$ labeled by column index $I$ and channel index $\alpha$. Such solutions will be combined to form the log-derivative matrix ${\mathbf{Y}}_{a}={\mathbf{F}}_{a}^{\prime}{\mathbf{F}}_{a}^{-1}$.

Our main equation (22) specialized to an element with $P_{m}$ points can now be rearranged as an algorithm expressing the value of ${\mathbf{Y}}_{P_{m}}$ on the right-end of the element to a known input value ${\mathbf{Y}}_{1}$ assigned on the left-end point. This task can be accomplished by right multiplications by ${\mathbf{F}}_{P_{m}}^{-1}$ to give after simple algebra :

[TABLE]

With $\bar{\mathbf{F}}_{a}\equiv{\mathbf{F}}_{a}{\mathbf{F}}_{P_{m}}^{-1}$ determined at first $(P_{m}-1)$ points, the remaining equation at last point

[TABLE]

determines the final log-derivative through a series of sparse matrix multiplications. The log-derivative ${\mathbf{Y}}_{P_{m}}$ can then be used as entry for the calculation in next element.

II.4 Error control

A major advantage of the spectral element method is the possibility to estimate precisely the numerical error by a posteriori analysis of the calculated solution. To this aim, let us consider the solution wavefunction restricted to the $m$-th sector and suppose for the sake of simplicity that the angular basis of dimension $N_{m}$ is constant within the sector. The discretized solution at the $P_{m}$ points in the sector for the different channels is therefore represented by $P_{m}\times N_{m}$ elements noted $F_{p\alpha}$.

The error is estimated by first performing an orthogonal transformation from the grid basis to the polynomial basis of Legendre polynomials $P_{n}(R)$ defined in the $[R_{1}^{(m)},R_{P_{m}}^{(m)}]$ interval through a coordinate transformation in the same fashion as in Eq. (3). Assuming the polynomials normalized, the transformation matrix reads explicitly $O_{np}=P_{n}(R_{p})/\sqrt{w_{p}^{(m)}}$. For each channel $\alpha$ the transformation $\mathbf{O}$ gives the set of pseudospectral coefficients ${\tilde{F}}_{n\alpha}$ of the solution expanded on the Legendre basis as ${\tilde{F}}_{n\alpha}=\sum_{n=1}^{P_{m}}O_{np}F_{p\alpha}$. The main point is that the convergence of the Legendre polynomial series is superalgebraic, at least for sufficiently regular solutions Boyd (2001). The size of last calculated coefficients ${\tilde{F}}_{P_{m}\alpha}$ is therefore a reliable estimate of the remainder of the series, i.e. of the numerical truncation error in each channel. If the error is larger (smaller) that a given tolerance criterion one can either reduce (increase) the element size or increase (reduce) the polynomial order $P_{m}$. As recognized at the birth of the so-called $hp$-methods, the optimal strategy to guarantee an exponential accuracy of the calculated solution consists in increasing $P_{m}$ in the regions where the latter is regular and in decreasing the element size in the regions where it is irregular Babuška and Guo (1992).

We will show in the next section a series of numerical experiments for both scattering and bound state calculations. We limit ourselves to a relatively simple ro-vibrational model with a purely diabatic basis in order to make the numerical convergence analysis as plane as possible. Since as most usual in molecular physics the solution is regular we fix the same polynomial order in all elements and study the behavior of selected observables as a function of both the element size and polynomial order.

III Numerical tests

We perform numerical tests of efficiency and accuracy of the algorithm on the Rb2He trimer, a system for which bound states and ultra-cold scattering properties have already been studied in our group Guillon et al. (2012); Viel and Launay (2014).

For the description of Rb2He, the $\vec{R}$ and $\vec{r}$ Jacobi vectors are used. The corresponding Hamiltonian in the space fixed frame reads Arthurs and Dalgarno (1960)

[TABLE]

with $R$ and $r$ the Rb2 – He and Rb2 distances, \mu_{\text{Rb{}{2}--He}} and \mu_{\text{Rb{}{2}}} the associated reduced masses, ${\vec{L}}$ and ${\vec{j}}$ the angular momenta and $\hat{V}$ the potential term taken from Ref. Viel and Launay, 2014 limited to the 2-body part. The generic $\Omega$ coordinates introduced in Sec. II correspond for this system to five spatial coordinates, namely, $r$, $\hat{r}$ and $\hat{R}$ that reduce to three when fixing the total angular momentum quantum numbers $J$ and $M$. The basis functions used to represent $\Phi_{a}(\Omega)$ in Eq. (21) are taken identical for all sectors. Equation (21) reads for this specific case

[TABLE]

where $\chi_{vj}(r)$ are the rovibrational eigenstates of the Rb2 diatomic and $Y_{j\ell}^{JM}$ the coupled spherical harmonics Arthurs and Dalgarno (1960).

III.1 Scattering states

For the calculations, we vary $R$ from 4 to 120 $a_{0}$, we use Rb${}_{2}(v=1,j=0)$ as the initial state for the collision, and we impose $R$-matrix boundary conditions. The linear system Eq. (25) is solved using the PARDISO package Petra et al. (2014a, b) included in the MKL library MKL . This state-of-the art direct solver determines the solution of a sparse linear system in (number of nonzero elements)3/2 operations. In our tests, we find that total memory used by PARDISO is about five times larger than the memory required to store the nonzero elements of the discretized Hamiltonian. The $K$-matrix, extracted from the matching procedure in Eq. (29) performed at $R=120$ $a_{0}$, is diagonalized to compute the eigenphasesum

[TABLE]

with $\eta_{i}$ being the $n_{op}$ eigenvalues of the open-open part of the $K$-matrix. Note that while the eigenphasesum is a function of collision energy and depends on the partial wave $J$ considered, the corresponding indices have been dropped for ease of notation.

For the accuracy tests, we used a collisional energy of 1 K above the $v=1,j=0$ initial state and focussed on the $J=1$ partial wave. Basis functions with up to $v=4,j=24$ quantum numbers are included, resulting in 125 channels, 32 of which are energetically open at the considered collision energy. A fixed Lobatto order $P$ is used for all $M$ elements used for the discretization of the $[R_{1},R_{L}]=[4,120]$ (in $a_{0}$) interval. All elements are taken of the same length, noted $h$ hereafter, with

[TABLE]

Figure 3 presents the variation of $\delta$ as a function of $h$ in log-log scale when $h$ is systematically divided by two.

Assuming a dependence of $\delta$ of the form

[TABLE]

the power $\alpha$ is directly accessible by the slope of the variation of the quantity $\Delta\delta=\delta(h)-\delta(h/2)$ as a function of $h$ in log-log scale, even when the exact value of $\delta_{0}$ is unknown. A comparison with the resolution of the coupled equations performed by the Johnson log-derivative propagator Johnson (1973) is also presented in the figure. The results for the Lobatto order $P=3$ case are strictly identical to the ones obtained with the Johnson propagator method. One can indeed show formally that solution of Eq. (36) with two Lobatto points followed by application of (37) gives exactly the same result as the three-points Johnson recursion. However, the two algorithms should not be viewed as equivalent, in the sense that the first half propagation step in Johnson’s method is not equivalent to solving our equation (36) for two points.

For these cases, the known $\alpha=4$ value can be read from the linear curves in the figure. For $\log_{10}(h)<-2.2$ the accuracy on $\delta$ can no longer be improved by a reduction of $h$ when using the Johnson propagator. For $P=3$ Lobatto case, the memory requirement prohibits the computation at such small $h$. For the $P=4$ to $P=7$ Lobatto orders, a slope $\alpha=2P-2$ is obtained. When increasing $P$, lower values of the absolute error are obtained for specific values of $h$. For example a $10^{-10}$ accuracy is reached for $P=7$ and $h\sim 0.4$ $a_{0}$. We stress that $\delta$ is a very sensitive quantity and such absolute error value on $\delta$ corresponds to the usually observed accuracy on rate coefficients calculations, much easier to converge.

An analysis of the CPU time needed for given accuracies is presented in Tab. 1.

The table presents as a function of the Lobatto order $P$ the CPU in second needed to reach a given accuracy of $\log_{10}(\Delta\delta)$. The corresponding $h$ are also listed in the table. For these calculations we increase $J$ to 10 for which the number of channels increases to 565 channels with 129 of them energetically open. The CPU given in the table corresponds to the resolution of the equations for one value of the collisional energy after an initialization step which is energy independent and thus to be performed only once if multiple collisions energies are considered. The table shows that for a given accuracy, one gains in increasing the Lobatto order at least to the tested orders. The comparison of the CPU time needed by the Johnson Johnson (1973) log-derivative propagator is clearly in favor of this last approach when low accuracy is required. However the situation changes when high accuracy is needed. In the present test, a $P=7$ computation is always more efficient than the Johnson version even for a single collision energy. When using the Gauss-Lobatto discretization, the CPU time requirement of the collisional energy dependent step is closely related to the number of integration points $L$ as underlined by Fig. 4.

A roughly linear dependence of the CPU time as a function of the number of integration points is found for the Lobatto orders we tested. At a given number of integration points, the general trend is an increase of the CPU time with the Lobatto order. Some exceptions are found like the $P=6$ case presented in Fig. 4 which turns out to be cheapest calculation with respect to CPU time for all numbers of integration points from 300 to 1000. We infer that this is due to particularities in the sparseness structure of the matrices handled by PARDISO. Similar but less marked exceptions have been found for the $J=10$ and $J=14$ partial wave computations.

III.2 Bound states

With the appropriate boundary conditions built in Eq. (24), bound states of the triatomic Rb2He can be determined using the same discretized Hamiltonian. We solve the sparse eigenvalue problem Eq. (24) using the density-matrix-based algorithm FEAST, a package included in the MKL library based on a contour representation of the resolvent in the complex plane Polizzi (2009). Internally, FEAST solves a series of sparse linear systems using a user-defined subroutine, PARDISO in our case.

We focus on the $J=2,\Pi=-$ partial wave for which a single bound state is found. This state is weakly bound with respect to the Rb2 + He asymptote and an enlarged $R$ box with $R\in[4,236]$ $a_{0}$ is used. The converged energy is -12.63 milliK below the Rb2 + He asymptote. Fig. 5 presents the evolution of the relative error on the computed energy when increasing by two the element size $h$ for various Lobatto order $P$ in log-log scale.

For the lowest $P$ orders presented, the $h^{2P-2}$ behavior is retrieved. For $P=5$ numerical noise increases when reducing $h$. This is partly due to the FEAST algorithm which implies an iterative procedure with two kinds of internal convergence criteria. One criterium tests the evolution of the energies from one iteration to the next one and the second one is a maximum number of iterations. The data presented have been obtained with a $10^{-14}$ value and 50 iterations maximum for these two FEAST parameters.

IV Conclusions

In summary, we have explored the numerical performance of the spectral-element method in multichannel quantum dynamics. Combination of the spectral-element discretization with purely diabatic or diabatic-by-sector bases leads to a highly sparse representation of the Hamiltonian. This results in significant memory saving for the bound state problem as compared for instance to the scaled DVR approach Tiesinga et al. (1998); Kokoouline et al. (1999). Regarding the scattering problem, accuracy is significantly less limited by round-off errors in the spectral element approach than in popular propagation methods and the corresponding computation time is advantageous when the required accuracy is high.

In perspective, it may be interesting to test iterative rather than direct algorithms to solve the scattering linear system for the discretized Schrödinger equation, in particular when boundary conditions of Eq. (34) are imposed to obtain a single column of the scattering matrix. In this case, if iterative solvers turned out to perform better than (number of nonzero elements)3/2 one might be able to overcome the (number of channel)3 unfavorable computational cost scaling presented by time-independent calculations as compared to time-dependent calculations.

Acknowledgements.

This work is supported by the Agence Nationale de la Recherche (Contract COLORI No. ANR-12-BS04-0020-01).