Domain decomposition method for the N-body time-independent and time-dependent Schrödinger equations

Abstract

This paper is devoted to the derivation of a pleasingly parallel Galerkin method for the time-independent N-body Schrödinger equation, and its time-dependent version modeling molecules subject to an external electric field (Bandrauk 1994; Bandrauk et al., J. Phys. B-Atom. Mol. Opt. Phys. 46(15), 153001, 2013; Cohen-Tannoudji et al. 1992). In this goal, we develop a Schwarz waveform relaxation (SWR) domain decomposition method (DDM) for the N-body Schrödinger equation. In order to optimize the efficiency and accuracy of the overall algorithm, (i) we use mollifiers to regularize the singular potentials and to approximate the Schrödinger Hamiltonian, (ii) we select appropriate orbitals, and (iii) we carefully derive and approximate the SWR transmission conditions. Some low-dimensional numerical experiments are presented to illustrate the methodology.

Appendices

Appendix A: Local orbital construction

In this appendix, we detail the construction of the local orbitals. We consider the two cases, for a given subdomain $D_{i}$, with $2\leq i \leq L-1$.

• $D_{i}$contains a nucleus singularity. Few subdomains belong to this first category, in particular when we are interested in the time-dependent Schrödinger equation for intense field-particle interaction. In that case, mollifiers will allow for an arbitrarily accurate smoothing of the nucleus singularities. Notice that in 1-d, the singularity treatment is different than in 3-d. Indeed in the latter case, we benefit from the fact that a Coulomb potential, up to a multiplicative constant is a fundamental solution to Poisson’s equation. This property allows for an accurate and efficient treatment of the localized orbitals. The Coulomb potential is then approximated by a smooth function $G_{\epsilon }$, thanks to mollifiers $B_{\epsilon }$ as defined in [1], and such that:

  $$\begin{array}{@{}rcl@{}} G_{\epsilon}(\mathbf{x}) = \frac{1}{4\pi}\left( V*B_{\epsilon}\right)(\mathbf{x}) \end{array} $$

  (A.1)

  where $V(\mathbf {x}) = -1/|\mathbf {x}|$, which also satisfies

  $$\begin{array}{@{}rcl@{}} 4\pi \triangle G_{\epsilon}(\mathbf{x}) = B_{\epsilon}(\mathbf{x}). \end{array} $$

  (A.2)

  As a consequence, a smooth approximation of the Coulomb potential V using (A.2), can be constructed with $G_{\epsilon }\rightarrow _{\epsilon \rightarrow 0}V$ in $\mathcal {D}^{\prime }({\mathbb {R}}^{3})$. Notice that this property is also fundamental for efficiently computing the 6-dimensional integrals in order to construct the global discrete Hamiltonian [1]. In 1-d, the fundamental solution of the Poisson equation is $|x|$ and the property just mentioned does not occur anymore. Instead, we directly compute G_ε using (A.1) with $B_{\epsilon }$ defined by:

  $$\begin{array}{@{}rcl@{}} B_{\epsilon}(x)=\left\{ \begin{array}{cc} \frac{1}{\epsilon}\sigma_{M(1)}\left( 1-\left( \frac{x}{\epsilon}\right)^{2}\right)^{M}, & |x|\leq \epsilon,\\ 0, & |x|>\epsilon \end{array} \right. \end{array} $$

  (A.3)

  where M refers to the order of the mollifier and the scaling factors $\sigma _{M}(1)$ are explicitly defined in [1]. For instance, for $\sigma (1)= 3/4$, $\sigma (2)= 15/16$ we represent in Fig. 18 for a unique domain $B_{\epsilon = 0.5}$ (Left) and $G_{\epsilon = 0.5}$ (Right). In particular it is proven in [1], that for any smooth function f

  $$\begin{array}{@{}rcl@{}} \|f-f*B_{\epsilon}\|_{2} = \sum\limits_{k \geq 1}c_{k}\epsilon^{2k} \end{array} $$

  for some positive sequence $\{c_{k}\}_{k}$. Once $G_{\epsilon }$ is computed, we introduce a barrier potential as in [3]

  $$\begin{array}{@{}rcl@{}} V_{\text{b}}(x) = s_{\epsilon_{b}}(x-x_{b})V_{\infty} \end{array} $$

  (A.4)

  where i) the smooth function $s_{\epsilon _{b}}$ is equal to 0 for $x<x_{b}-\epsilon _{b}/2$ and 1 for $x>x_{b}-\epsilon _{b}/2$, ii) $\epsilon _{b}>0$, and iii) $V_{\infty }$ and $x_{b}$ the center of the barrier potential, are imposed. The support of the localized orbitals is then $(x_{c_{i}}-x_{b},x_{c_{i}}+x_{b})$, where $x_{c_{i}}$ denotes the coordinates of the center of the subdomain $D_{i}$. We typically choose $x_{b} > |D_{i}|/2$ to ensure that the localized orbitals are not null at $D_{i}$’s boundary. A contrario, taking $x_{b}$ too large will lead to a loss of computational efficiency due to a large localized orbital support. In $D_{i}$, we then solve the following one-dimensional one-electron eigenvalue problem

  $$\begin{array}{@{}rcl@{}} \left( -\frac{1}{2}\partial_{x}\left( a_{\epsilon}(x)\partial_{x}\right) + G_{\epsilon}(x-\overline{x}_{A}) + G_{\epsilon}(x-\overline{x}_{B}) + V_{\text{b}}(x-x_{c_{i}})\right){\varphi_{l}^{i}}(x) = \lambda^{i}_{l,\epsilon}{\varphi_{l}^{i}}(x) \end{array} $$

  where $1\leq l\leq M_{i}$ (resp. $2 \leq i\leq L-1$) is the orbital (resp. subdomain) index and where $a_{\epsilon }(x):= 1-s_{\epsilon _{b}}(x-x_{b})$. Notice that the choice of the localized orbitals is motivated by physical considerations. When we are interested in field-particle interaction, for subdomains containing the nuclei, we will select the localized orbitals corresponding to the lower energy states, as they will be predominant in the overall wavefunction in the vicinity of the nucleus singularities.

  Fig. 18

  

  (Left) Mollifiers B_0.5 (M = 1,2), (Right) and G_0.5 for d = 1

• $D_{i}$does not contain any nucleus singularity. In that case, the regularization of the Coulomb potential through mollifiers is naturally useless. The localized orbitals are then directly obtained by solving

  $$\begin{array}{@{}rcl@{}} \left( -\frac{1}{2}\partial_{x}\left( a_{\epsilon}(x)\partial_{x}\right) - \frac{1}{|x-\overline{x}_{A}|} - \frac{1}{|x-\overline{x}_{B}|}+ V_{\text{b}}(x)\right){\varphi_{l}^{i}}(x) = \lambda^{i}_{l,\epsilon}{\varphi_{l}^{i}}(x) \end{array} $$

  Similarly to the previous case (subdomain containing the nuclei), the selected localized orbitals will strongly depend on the relative position of the nuclei $/ D_{i}$. Alternatively, for those subdomains, we can use local Gaussian basis functions, as described in Section 2.1.

Once the localized orbitals are computed, we can construct the discrete Schrödinger Hamiltonian. The approach which is proposed below will benefit from i) the fact that the localized orbitals are selected accordingly to the position of the nuclei, ii) the compact support of the LO’s and iii) their orthogonality property (more specifically their extension by 0 to all ${\Omega }$). This last point necessitates some precisions. First, we notice that by construction for any $i \in \{2,\cdots ,L-1\}$, the supports of $\left \{{\varphi _{l}^{i}}\right \}_{1\leq l\leq M_{i}}$ and of $\big \{{\varphi _{m}^{j}}\big \}_{1\leq m\leq M_{j}}$ with $j \neq i-1,i,i + 1$ are disjoint, so that these LO’s are trivially orthogonal. By construction, the LO’s $\left \{{\varphi _{l}^{i}}\right \}_{l}$ of any $D_{i}$ are also orthogonal to each other. For $j=i-1$ or $j=i + 1$, the orthogonality of the LO’s $\left \{{\varphi _{l}^{i}}\right \}_{1\leq l\leq M_{i}}$ and of $\big \{{\varphi _{m}^{j}}\big \}_{1\leq m\leq M_{j}}$ is not, a priori, ensured. However, by construction for any $1\leq l \leq M_{i}$ and $1\leq m \leq M_{i\pm 1}$

$$\begin{array}{@{}rcl@{}} \Big|\text{Supp}\left( {\varphi_{l}^{i}}\cap\varphi^{i\pm 1}_{m}\right)\Big| \leq 2x_{b}-\big|D_{i} \cup D_{i\pm 1}\big| . \end{array} $$

Then, as the LO’s (smoothly) vanish at the boundary of their support, for $x_{b}-|D_{i,i\pm 1}|/2$ small enough, we expect that ${\int }_{{\mathbb {R}}}{\varphi _{l}^{i}}(x)\varphi ^{i\pm 1}_{m}(x)dx$ to be small. For $L = 2$ and $x_{b}= 8$ (which is relatively very large) and (L =) 5 subdomains, we represent for subdomain $D_{i = 0}=(-10,10)$ (resp. $D_{i = 1}=(-2,18)$) ${\varphi _{l}^{0}}$ (resp. ${\varphi ^{1}_{l}}$), for $l = 1,\cdots ,4$ in Fig. 19.

Fig. 19

First 4 eigenstates LO’s $\phi _{l,\epsilon }^{3}$ (Left), $\phi _{l,\epsilon }^{4}$ (Right) for l = 1,⋯ ,4

Appendix B: Antisymmetric wavefunction reconstruction

In order to construct an antisymmetric wavefunction, we propose a specific decomposition of ${\mathbb {R}}^{dN}$ with antisymmetric local basis functions. In order to simplify the presentation, we will assume that i) the subdomains $\{{\Omega }_{i}\}_{i \in \{1,\cdots ,L^{dN}\}}\in {\mathbb {R}}^{dN}$ are hypercubes of identical size and $L\in 2{\mathbb {N}}+ 1$, and ii) there is no overlap between the subdomains. We first define:

Definition Appendix B.1

We denote by $\sigma (i;p,q) \in \{1,\cdots ,L^{dN}\}$ the subdomain index such that for (x₁,⋯ ,x_p,⋯ ,x_q,⋯ ,x_N) ∈Ω_i: $(\mathbf {x}_{1},\cdots ,\mathbf {x}_{q},$$\cdots ,\mathbf {x}_{p},\cdots ,\mathbf {x}_{N}) \in {\Omega }_{\sigma (i;p,q)}$. The index i refers to the subdomain index, and (p,q) to the permutation coordinate indices.Notice that we naturally have $\sigma \left (\sigma (i;p,q);q,p\right )=i$ and that $\sigma \left (\sigma (i;p,q);q,p\right )$ is unique as there is no subdomain overlap. Antisymmetry of the wavefunction would then occur if for all $i \in \{1,\cdots ,L^{dN}\}$, the local basis $\left \{{v^{i}_{l}}\right \}_{1\leq l\leq K_{i}}$ coincides with the local bases $\left \{v^{\sigma (i;p,q)}_{l}\right \}_{1\leq l \leq K_{\sigma (i;p,q)}}$, for all $(p,q) \in \{1,\cdots ,N\}^{2}$. We define

$$\begin{array}{@{}rcl@{}} {\Sigma}(i)=\left\{\sigma(i;p,q)\in \{1,\cdots,L^{dN}\}, \forall (p,q)\in \{1,\cdots,N\}^{2}\right\} \end{array} $$

As there is no overlap between subdomains, ${\Sigma }(i)$ is actually a singleton. In order to guarantee the antisymmetry of the overall wavefunction, at any time and any Schwarz iteration, we proceed as follows.

Antisymmetric wavefunction

For any i ∈{1,⋯ ,L^dN} and for any $l \in {\Sigma }(i)$:

• we solve a (local) TDSE on ${\Omega }_{i}$ in the form

  $$\begin{array}{@{}rcl@{}} \psi^{(k)}_{i}(\cdot,t) = \sum\limits_{l = 1}^{K_{i}}c_{l}^{i,(k)}(t){v^{i}_{l}} \end{array} $$

  (B.1)

• and we deduce:

  $$\begin{array}{@{}rcl@{}} \psi^{(k)}_{{\Sigma}(i)}(\cdot,t) = -\sum\limits_{l = 1}^{K_{{\Sigma}(i)}}c_{l}^{i,(k)}(t)v^{{\Sigma}(i)}_{l} \end{array} $$

  (B.2)

The global wavefunction is then reconstructed according to the algorithm derived in Section 3.3 without overlap. We deduce by construction, the following proposition.

Proposition Appendix B.1

The reconstructed wavefunction $\psi ^{(k)}$with (B.1) and (B.2) is antisymmetric.

Proof

The proof is trivial. For any $(\mathbf {x}_{1},\cdots ,\mathbf {x}_{p},\cdots ,\mathbf {x}_{q},\cdots ,\mathbf {x}_{N}) \in {\Omega }_{i}$, we have (x₁,⋯ ,x_q,⋯ ,x_p,⋯ ,x_N) ∈Ω_Σ(i), then

$$\begin{array}{@{}rcl@{}} \psi^{(k)}(\mathbf{x}_{1},\cdots,\mathbf{x}_{p}, \cdots,\mathbf{x}_{q},\cdots, \mathbf{x}_{N},t) & =& \psi_{i}^{(k)}(\mathbf{x}_{1},\cdots,\mathbf{x}_{p}, \cdots,\mathbf{x}_{q},\cdots, \mathbf{x}_{N},t)\\ \\ & = & -\psi_{{\Sigma}(i)}^{(k)}(\mathbf{x}_{1},\cdots,\mathbf{x}_{q}, \cdots,\mathbf{x}_{p},\cdots, \mathbf{x}_{N},t)\\ \\ & =& -\psi^{(k)}(\mathbf{x}_{1},\cdots,\mathbf{x}_{q}, \cdots,\mathbf{x}_{p},\cdots, \mathbf{x}_{N},t). \end{array} $$

□

Discrete local Hamiltonian construction with local Slater’s determinants

The construction of discrete local Hamiltonians $\widetilde {\textbf {H}}_{i}$ is similar to the procedure described in Section 2.2 and Appendix A, which is an application of [2], in a DDM framework. Due to the antisymmetry constraint detailed above, we however need additional specifications. The strategy presented in Section 2.2, and extended to d dimensions and N particles allows for the construction for any subdomain Ω_i, $1\leq i \leq L^{dN}$, of the SLO’s $\left \{{\varphi _{l}^{i}}(\mathbf {x})\right \}_{l}$. These SLO’s are smooth, have compact support and possess orthogonality properties, which are described at the end of Appendix A. We here summarize the explicit construction of a local Hamiltonian $\widetilde {\textbf {H}}_{i}$, $1 \leq i \leq L^{dN}$, say for Ω_i. Following the notations used above, we need to compute:

$$\begin{array}{@{}rcl@{}} I_{mpqr}^{ijkl} &= &{\int}_{{\mathbb{R}}^{dN}} {\varphi^{i}_{m}}(\mathbf{x}){\varphi^{j}_{p}}(\mathbf{x})\frac{{\varphi^{k}_{q}}(\textbf{s}){\varphi^{l}_{r}}(\textbf{s})}{|\textbf{s}-\mathbf{x}|}d\mathbf{x}d\textbf{s} \\ &= &{\int}_{\text{Supp}{\varphi_{m}^{i}}\cap \text{Supp}{\varphi_{p}^{j}}\cap \text{Supp}{\varphi_{q}^{k}}\cap \text{Supp}{\varphi_{r}^{l}}}{\varphi^{m}_{j}}(\mathbf{x}){\varphi^{j}_{p}}(\mathbf{x}){\Phi}_{qr}^{kl}(\mathbf{x})d\mathbf{x}. \end{array} $$

where

$$\begin{array}{@{}rcl@{}} {\Phi}_{qr}^{kl}(\mathbf{x}) = {\int}_{\text{Supp}{\varphi^{k}_{q}}\cap \text{Supp}{\varphi^{l}_{r}}}\frac{{\varphi_{q}^{k}}(\textbf{s}){\varphi_{r}^{l}}(\textbf{s})}{|\textbf{s}-\mathbf{x}|}d\textbf{s} \end{array} $$

can be achieved using mollifiers:

$$\begin{array}{@{}rcl@{}} {\Phi}_{qr}^{kl}(\mathbf{x}) \approx {\int}_{\text{Supp}{\varphi^{k}_{q}}\cap \text{Supp}{\varphi^{l}_{r}}}{\varphi_{q}^{k}}(\textbf{s}){\varphi_{r}^{l}}(\textbf{s})B_{\epsilon}(\textbf{s}-\mathbf{x})d\textbf{s} \end{array} $$

or alternatively, if $d = 3$

$$\begin{array}{@{}rcl@{}} -\triangle {\Phi}_{qr}^{kl} = 4{\pi\varphi_{q}^{k}}{\varphi_{r}^{l}} \end{array} $$

Then, for $m=l,l\pm 1$

$$\begin{array}{@{}rcl@{}} I_{lm}^{ij} = {\int}_{\text{Supp}{\varphi^{i}_{l}}\cap \text{Supp}{\varphi_{m}^{j}}}\frac{1}{2}{\varphi_{l}^{i}}(\mathbf{x}){\triangle\varphi_{m}^{j}}(\mathbf{x})-\sum\limits_{A = 1}^{P}\frac{Z_{A}}{|\mathbf{x}-\mathbf{x}_{A}|}{\varphi_{l}^{i}}(\mathbf{x}){\varphi_{m}^{j}}(\mathbf{x})d\mathbf{x}. \end{array} $$

Notice that for constructing $\textbf {Q}_{i}^{x,y,z}$, we need to compute

$$\begin{array}{@{}rcl@{}} J^{ij}_{lm} = {\int}_{\text{Supp}{\varphi^{i}_{l}}\cap \text{Supp}{\varphi_{m}^{j}}}\mathbf{x}{\varphi_{l}^{i}}(\mathbf{x}){\varphi_{m}^{j}}(\mathbf{x})d\mathbf{x} \end{array} $$

which does not present any additional difficulty compared to one-domain problems, see again [3]. Efficiently construction of the local Hamiltonians can be performed using the strategy presented in [2]. Once the matrix local Hamiltonians are constructed, we can determine the LSD’s and solve the time-independent and -dependent Schrödinger equations.

Appendix C: Computational complexity analysis of the SWR-DDM solver for N-body equation

In this appendix, we analyze the computational complexity of the SWR method for both the time-independent and time-dependent cases. An obvious consequence of the use of a SWR/FCI method with orbital basis functions, is that the number of degrees of freedom (dof) is expected to be much smaller compared to finite difference/volume methods (FDM/FVM) or low degree finite element methods (FEM). For instance, in a subdomain ${\Omega }_{i}$, a cell center FVM (Q₀) consists of choosing the basis functions as $\left \{\textbf {1}_{{V^{i}_{l}}}(\mathbf {x}_{1},\cdots ,\mathbf {x}_{N})/|{V^{i}_{l}}|\right \}_{1\leq l\leq N_{i}}$ with finite volumes ${V^{i}_{l}}$, such that $\cup _{l = 1}^{N_{i}}{V_{l}^{i}}=\tau _{h}({\Omega }_{i})$. That is:

$$\begin{array}{@{}rcl@{}} \psi_{i}^{(k)}(\mathbf{x}_{1},\cdots,\mathbf{x}_{N},t) = \sum\limits_{l = 1}^{N_{i}}\frac{1}{|{V^{i}_{l}}|}\textbf{1}_{{V_{l}^{i}}}{c^{i}_{l}}(t), (\mathbf{x}_{1},\cdots,\mathbf{x}_{N}) \in {\Omega}_{i}. \end{array} $$

Although $Q_{0}$-basis functions are very simple, and that the corresponding matrices are sparse, in order to get a precise description of the wavefunction a “very” large number $N_{i}$ on ${\Omega }_{i}$, of finite volumes is necessary. Slater’s orbitals $\left \{{v_{l}^{i}}\right \}_{1\leq l\leq K_{i}}$, defined above, would typically contain a very large number of finite volumes ${V_{l}^{i}}$ (N_i ≫ K_i). The consequence is that, although FVM-linear systems are much sparser than Galerkin-FCI systems, they are also of much higher dimension.

Below, we study the overall computational complexity (CC) and the scalability of the SWR-DDM in d-dimension and for N particles. The analysis will be provided for both the stationary and unstationary N-body Schrödinger equations.

Computational complexity and scalability for the time-independent N-body equation

From now on, we assume that the computational domain is decomposed in L^dN-subdomains, $\{{\Omega }_{p}\}_{1\leq p \leq L^{dN}}$, and that $K_{p}$ (resp. $\mathcal {K}_{p}$) with $1 \leq p \leq L^{dN}$, local basis functions are selected (resp. the number of degrees of freedom) per subdomain. We denote by $K_{\text {Tot}}:={\sum }_{p = 1}^{dN} K_{p}$ (resp. $\mathcal {K}_{\text {Tot}}:={\sum }_{p = 1}^{dN} \mathcal {K}_{p}$) the total number of local basis functions (resp. degrees of freedom). In the following, we will assume for simplicity that there is a fixed number of local basis functions per subdomain, that is $\mathcal {K}_{p} \approx \mathcal {K}_{\text {Tot}}/L^{dN}$, for any $1 \leq p \leq L^{dN}$. Instead of dealing with a full discrete Hamiltonian in $M_{\mathcal {K}_{\text {Tot}}}({\mathbb {R}})$, we then rather deal with $L^{dN}$ local discrete Hamiltonians in $M_{\mathcal {K}_{\text {Tot}}/L^{dN}}({\mathbb {R}})$. The very first step then consists of constructing the local basis functions, then of the $L^{dN}$ local Hamiltonians, which is pleasingly parallel (perfect distribution of the integral computations). We focus on the complexity and scalability for computing the eigenenergies from these (local or global) discrete Hamiltonians.

We then recall the main ingredients necessary to study the computational complexity analysis of the SWR method presented in Section 3 for solving the time-independent Schrödinger equation using the NGF method. We have decomposed the spatial domain ${\Omega } \subset {\mathbb {R}}^{dN}$ in $L^{dN}$ subdomains and solve an imaginary-time-dependent Schrödinger equation (or real-time normalized heat equation) on each subdomain. At a given Schwarz iteration $k \in {\mathbb {N}}$, we denote by $T^{(k)}_{p}$ (resp. $n_{p}^{(k)}$) the imaginary convergence time (resp. number of time iterations to converge) for the NGF method in the subdomain ${\Omega }_{p}$, where $1\leq p \leq L^{dN}$. In addition, each imaginary time iteration requires $\mathcal {O}(\mathcal {K}_{p}^{\beta ^{\text {(S)}}_{p}})$ operations, where $1 < \beta ^{\text {(S)}}_{p} < 3$ (coming from sparse linear system solver). The index $(\text {S})$ refers to the stationary Schrödinger equation. We denote by $k^{\text {(cvg)}}$, the total number of Schwarz iterations to reach convergence, as described in Section 3. Notice that $k^{\text {(cvg)}}$ is strongly dependent on the type of transmission conditions [8]. We get

Proposition Appendix C.1

The computational complexity CC$^{\text {\text {(S)}}}_{\text {SWR}}$ofthe overall SWR-DDM method described in Section 3.2 for solving the Schrödingerequation in the stationary case is given by

$$\begin{array}{@{}rcl@{}} \text{CC}^{\text{(S)}}_{\text{SWR}} = \mathcal{O}\left( \sum\limits_{k = 1}^{k^{\text{(cvg)}}}\sum\limits_{p = 1}^{L^{dN}}n_{p}^{(k)}\mathcal{K}_{p}^{\beta^{\text{(S)}}_{p}}\right). \end{array} $$

(C.1)

Assuming that $\beta ^{\text {(S)}}_{p}$,(resp. $n_{p}^{(k)}$)isps-independent (that is subdomain-independent), and thendenoted $\beta ^{\text {(S)}}$(resp. $N^{(k)}$),we have

$$\begin{array}{@{}rcl@{}} \text{CC}^{\text{(S)}}_{\text{SWR}} = \mathcal{O}\left( \frac{\mathcal{K}_{\text{Tot}}^{\beta^{\text{(S)}}}}{L^{dN(\beta^{\text{(S)}}-1)}}\sum\limits_{k = 1}^{k^{\text{(cvg)}}} N^{(k)}\right). \end{array} $$

Thus

• assuming that the algorithm is implemented on a P-core machine, each core will deal with $\approx L^{dN}/P$ subdomains. The message passing load is dependent on the type of transmission conditions, but typically for classical or Robin SWR the communication will be very light. As a consequence an efficiency (T/PT_P) close to 1 is expected.

• we are dealing with $L^{dN}$ linear systems with approximately $\mathcal {K}_{\text {Tot}}/L^{dN}$ degrees of freedom (instead of a unique large system of $\mathcal {K}_{\text {Tot}}$ degrees of freedom if a huge discrete Hamiltonian was considered). As $\beta ^{\text {(S)}}$ is strictly greater than 1, we benefit from a scaling effect.

The SWR-DDM approach will be relevant in the starionary case, if typically

$$\begin{array}{@{}rcl@{}} \sum\limits_{k = 1}^{k^{\text{(cvg)}}} N^{(k)} \ll L^{dN(\beta^{\text{(S)}}-1)} . \end{array} $$

In order to satisfy this condition i) an implicit solver for the heat equation will allow for a faster convergence (the bigger the time step, the smaller $N^{(k)}$) of the NGF method, appropriate transmission conditions will allow for a minimization of $k^{\text {(cvg)}}$.

Computational complexity and scalability for the time-dependent N-body equation

Thanks to the SWR approach and as in the stationary case, the computation of the time-dependent equation, does not involve a full discrete Hamiltonian in $M_{\mathcal {K}_{\text {Tot}}}({\mathbb {C}})$, but rather $L^{dN}$ discrete local Hamiltonians in $M_{\mathcal {K}_{\text {Tot}}/L^{dN}}({\mathbb {C}})$. Notice that if we use the same Gaussian basis functions in each subdomain, the local potential-free Hamiltonians are identical in each subdomain, and has to be performed only once. The contribution from the interaction potential and laser field in the local Hamitonians, are however subdomain-dependent, but are diagonal operators. We assume that the TDSE is computed from time 0 to time $T>0$. An implicit scheme (L²-norm preserving) is implemented, which necessitates the numerical solution at each time iteration of a sparse linear system. We denote by $n_{T}$ the total number of time iterations to reach T, which will be assumed to be the same for both methods. We deduce the following proposition.

Proposition Appendix C.2

The computational complexity $\text {CC}^{\text {(NS)}}_{\text {SWR}}$ofthe overall SWR-DDM method described in Section 3.1 for solving the Schrödingerequation in the time-dependent case is given by:

$$\begin{array}{@{}rcl@{}} \text{CC}^{\text{(NS)}}_{\text{SWR}} = \mathcal{O}\left( n_{T}k^{\text{(cvg)}}\sum\limits_{p = 1}^{L^{dN}}\mathcal{K}_{p}^{\beta^{\text{(NS)}}_{p}}\right) \end{array} $$

(C.2)

with $1< \beta _{p}^{\text {(NS)}} < 3$, where theindex $\text {(NS)}$refersto the unstationary case. Assuming thatthe $\beta ^{\text {(NS)}}_{p}$isp-independent (that is subdomain-independent), anddenoted $\beta ^{\text {(NS)}}$,we get

$$\begin{array}{@{}rcl@{}} \text{CC}^{\text{(NS)}}_{\text{SWR}} = \mathcal{O}\left( n_{T}k^{\text{(cvg)}}\frac{\mathcal{K}_{\text{Tot}}^{\beta^{\text{(NS)}}}}{L^{dN(\beta^{\text{(NS)}}-1)}}\right). \end{array} $$

The SWR-DDM will then be relevant in the unstationary case, if

$$\begin{array}{@{}rcl@{}} k^{\text{(cvg)}} \ll L^{dN(\beta^{\text{(NS)}}-1)} . \end{array} $$