<figure><img src="image_1.png" /></figure>

HAL
open science

# Frozen Gaussian Approximation based domain decomposition methods for the linear Schrödinger equation beyond the semi-classical regime

Emmanuel Lorin, Xu Yang, Xavier Antoine

► To cite this version:

Emmanuel Lorin, Xu Yang, Xavier Antoine. Frozen Gaussian Approximation based domain decomposition methods for the linear Schrödinger equation beyond the semi-classical regime. Journal of Computational Physics, 2016, 315, pp.221-237. 10.1016/j.jcp.2016.02.035. hal-01244430

HAL Id: hal-01244430

https://hal.science/hal-01244430

Submitted on 7 Feb 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Frozen Gaussian Approximation based domain decomposition methods for the linear Schrödinger equation beyond the semi-classical regime

E. Lorin<sup>a,b</sup>, X. Yang<sup>c</sup>, X. Antoine<sup>d</sup>

<sup>a</sup>Centre de Recherches Mathématiques, Université de Montréal, Montréal, Canada, H3T 1J4

<sup>b</sup>School of Mathematics and Statistics, Carleton University, Ottawa, Canada, K1S 5B6

<sup>c</sup>University of California at Santa Barbara, Department of Mathematics, USA

<sup>d</sup>Institut Elie Cartan de Lorraine, Université de Lorraine, Inria Nancy-Grand Est, F-54506 Vandoeuvre-lès-Nancy Cedex, France

## Abstract

The paper is devoted to develop efficient domain decomposition methods for the linear Schrödinger equation beyond the semiclassical regime, which does not carry a small enough rescaled Planck constant for asymptotic methods (e.g. geometric optics) to produce a good accuracy, but which is too computationally expensive if direct methods (e.g. finite difference) are applied. This belongs to the category of computing middle-frequency wave propagation, where neither asymptotic nor direct methods can be directly used with both efficiency and accuracy. Motivated by recent works of the authors on absorbing boundary conditions [X. Antoine et al, J. Comput. Phys., 277 (2014), 268-304] and [X. Yang and J. Zhang, SIAM J. Numer. Anal., 52 (2014), 808-831], we introduce Semiclassical Schwarz Waveform Relaxation methods (SSWR), which are seamless integrations of semiclassical approximation to Schwarz Waveform Relaxation methods. Two versions are proposed respectively based on Herman-Kluk propagation and geometric optics, and we prove the convergence and provide numerical evidence of efficiency and accuracy of these methods.

**Keywords:** Semiclassical approximation, absorbing boundary condition, domain decomposition, linear Schrödinger equation, pseudodifferential operators.

## 1. Introduction

We are interested in computing the Time-Dependent linear Schrödinger Equation (TDSE) in the middle-frequency regime:

$$i\epsilon\frac{\partial\psi^{\epsilon}}{\partial t} = -\frac{\epsilon^2}{2}\Delta\psi^{\epsilon}(t, \mathbf{x}) + V(\mathbf{x}, t)\psi^{\epsilon}(t, \mathbf{x}), \mathbf{x} \in \mathbb{R}^d, t \in (0, T], \quad (1.1)$$

with $\psi^{\epsilon}(0, \mathbf{x}) = \psi_0^{\epsilon}(\mathbf{x}) \in L^2(\mathbb{R}^d, \mathbb{C})$ as the initial wave function, $\epsilon > 0$ as the rescaled Planck constant, $V$ as the potential function and $i = \sqrt{-1}$. The middle-frequency regime is calibrated as $\epsilon_0 \le \epsilon \ll 1$ whereas $\epsilon \le \epsilon_0$ is the high-frequency regime.

Direct numerical simulation of (1.1) is prohibitively expensive due to the small parameter $\epsilon \ll 1$. For example a mesh size of $O(\epsilon)$ is required when using the time-splitting spectral method [15], and the mesh size (and the time step as well) becomes as small as $o(\epsilon)$, if finite difference methods are used [37, 38]. In the meantime, $\epsilon \gg \epsilon_0$ is not small enough for semiclassical methods (e.g geometric optics, Gaussian beam methods, and frozen Gaussian approximation) to provide good accuracy [20, 28, 29, 34]. This motivates us to derive, analyze and present preliminary numerical performance of Semiclassical Schwarz Waveform Relaxation (SSWR) Domain Decomposition Methods (DDM) for the computation of (1.1) in the middle-frequency regime.

Schwarz Waveform Relaxation (SWR) methods are iterative fixed-point methods aiming to reconstruct global solutions from locally computed ones. Important efforts have been put on the development of Optimized Schwarz Waveform Relaxation (OSWR) methods for linear and nonlinear Schrödinger equations,
diffusion-advection-reaction equations, and Maxwell's equations [11, 12, 17, 18, 21, 24]. The principle of OSWR is to derive transmission conditions from Transparent or high-order Absorbing Boundary Conditions (TBCs/ABCs). The latter allow for an almost null spurious reflection of waves at boundaries of a computational domain. Such conditions were developed for many kinds of wave equations. We refer here to some of analytical and numerical works, mainly based on pseudodifferential operators (although other robust techniques exist such as perfectly matched layers (PMLs) [16]): wave equations [19], diffusion equation [22, 23], linear [1, 5, 7, 8, 9, 33], nonlinear Schrödinger equations [3, 6, 10, 30, 40, 42], Maxwell [2], Klein-Gordon and Dirac equations [13]. From the DDM point of view, and most particularly for OSWR-type methods, the use of TBC-based transmission conditions leads to fast converging algorithms (in a few fixed-point iterations). Therefore, most of the works dedicated to the development and discretization of TBCs benefit to the design of efficient and robust OSWR methods. In addition, the precise formal analysis of the rate of convergence of OSWR is also possible through the use of pseudodifferential operator techniques; see [11, 24] for instance. Although OSWR methods are in principle very attractive, they may nevertheless suffer from computational issues. Indeed, TBCs for the TDSE are defined through nonlocal (fractional) pseudodifferential operators which are computationally nontrivial and expensive to approximate, most particularly if one wants to preserve accuracy and stability properties of the interior schemes. More simple SWR techniques exist, such as the Classical Schwarz Waveform Relaxation (CSWR) method based on Dirichlet boundary conditions. Even if these methods are simple, they suffer from a slow convergence (and sometimes do not even converge), as discussed in [11, 24]. More elaborated SWR methods such as the ones based on optimized Robin boundary conditions [24] offer a good compromise (relatively fast convergence and small computational complexity).

In this paper, we develop Semiclassical Schwarz Waveform Relaxation (SSWR) methods for (1.1) in the middle-frequency regime. The originality primarily comes from the use of semiclassical solutions to the Schrödinger equation, more specifically the Frozen Gaussian Approximation (FGA) which provides a very accurate but relatively “inexpensive” ansatz in the semiclassical regime. This method was originally developed by Herman-Kluk (HK) [25] and was later mathematically analyzed in [39]. More recently, the HK formalism was used and analyzed to derive fast numerical solvers in the semiclassical regime for different classes of partial differential equations: scalar wave equations [34, 32], linear hyperbolic systems of conservation laws [35, 36], and ABCs for the linear Schrödinger equation in the semiclassical regime [43].

When the exact solution to the Schrödinger equation is known at subdomain interfaces, it is possible to derive extremely efficient DDMs. Although in general it is not possible to determine such information, some accurate approximations can be provided through FGA. In addition, as we will show, the FGA is well-adapted to parallel computing. The main purpose of this paper is to derive a SWR-DDM by using the semiclassical approximation for the Schrödinger equation in the middle-frequency regime, where a pure FGA or Geometric Optics Approximation (GOA) would be more appropriate than using the solution to the Schrödinger equation in quantum regime. The general strategy can then be summarized as follows: the Schrödinger equation is solved by using a SWR-DDM with FGA (or GOA)-based transmission conditions. Analysis and numerical simulations will be presented to show the main features of these original methods.

The rest of the paper is organized as follows. In Section 2, we design a FGA-based DDM in the semiclassical (high-frequency) regime. An analysis is provided to show the scalability of this DDM. Section 3 is dedicated to analyze the analogies of TBCs between the formalism of FGA and OSWR, based on which, in Section 4, we derive Herman-Kluk Schwarz Waveform Relaxation methods for the Schrödinger equation in the middle frequency regime. Alternative Schwarz Waveform Relaxation methods based on geometric optics are proposed in Section 5. Numerical experiments are provided in Section 6 to show the efficiency of the proposed methods, and we conclude in Section 7.

## 2. FGA-based domain decomposition method in the high-frequency regime

As a warmup, we first study the domain decomposition method (DDM) for solving the Schrödinger equation (1.1) in the high-frequency regime ($\epsilon \le \epsilon_0$) based on frozen Gaussian approximation (FGA), and then generalize it to middle-frequency regime in later sections.

We first recall that solving (1.1) by DDM often requires a Schwarz algorithm, which can be either derived from Dirichlet (CSWR) or transparent boundary conditions (OSWR) [12]. More specifically, let us
introduce two open sets $\Omega_{\eta}^{\pm}$ such that $\mathbb{R}^d = \Omega_{\eta}^{+} \cup \Omega_{\eta}^{-}$, where $\Omega_{\eta}^{+} \cap \Omega_{\eta}^{-}$ is in general non-empty with $\eta$ as a (usually small) parameter characterizing the overlapping region size. For example, for $d = 1$, one can choose $\Omega_{\eta}^{-} := (-\infty, \eta/2)$ and $\Omega_{\eta}^{+} := (-\eta/2, +\infty)$.

We denote by $\phi^{\epsilon, \pm}$ the solution to the time-dependent Schrödinger equation in $\Omega_{\eta}^{\pm}$. Then for any Schwarz iteration $k \ge 1$, the equation in $\Omega_{\eta}^{\pm}$ reads

$$ \left\{ \begin{array}{ll} i\epsilon\partial_t\phi^{\epsilon,\pm,(k)} = H_\epsilon\phi^{\epsilon,\pm,(k)} := \left(-\frac{\epsilon^2}{2}\Delta + V(\mathbf{x},t)\right)\phi^{\epsilon,\pm,(k)} & \text{in } \Omega_{\eta}^{\pm} \times (0,T), \\ B^{\epsilon,\pm}\phi^{\epsilon,\pm,(k)} = B^{\epsilon,\pm}\phi^{\epsilon,\mp,(k-1)} & \text{on } \Gamma_{\eta}^{\pm} \times (0,T), \end{array} \right. \qquad (2.1) $$

where $\phi^{\epsilon, \pm, (0)}$ are two given functions and $\Gamma_{\eta}^{\pm} = \partial\Omega_{\eta}^{\pm}$. For CSWR, $B^{\epsilon, \pm}$ is the identity operator and, for OSWR, $B^{\epsilon, \pm}$ is a nonlocal Dirichlet-to-Neumann-like (DtN) pseudodifferential operator. We refer to [12, 21, 24] for further reading.

Here we propose a FGA-based DDM by decomposition of the phase space rather than the real space. We introduce a Fourier Integral Operator $\mathcal{I}^{\epsilon}(t; a)$ with $t \ge 0$ and $a$ as a $(t, \mathbf{q}, \mathbf{p})$-dependent amplitude function. The FGA reads [43], for a given initial condition $\psi_0^{\epsilon}$ and $\mathbf{x} \in \mathbb{R}^d$,

$$ \mathcal{I}^{\epsilon}(t; a)\psi_0^{\epsilon}(\mathbf{x}) = \phi_{\text{FGA}}^{\epsilon}(t, \mathbf{x}) = \frac{1}{(2\pi\epsilon)^{3d/2}} \int_{\mathbb{R}^{2d}} a(t, \mathbf{q}, \mathbf{p})w^{\epsilon}(\mathbf{q}, \mathbf{p})e^{\frac{i}{\epsilon}(S(t, \mathbf{q}, \mathbf{p}) + \mathbf{P}(t, \mathbf{q}, \mathbf{p}) \cdot (\mathbf{x} - \mathbf{Q}(t, \mathbf{q}, \mathbf{p})))} \\ \times e^{-\frac{1}{2\epsilon}|\mathbf{x} - \mathbf{Q}(t, \mathbf{q}, \mathbf{p})|^2} d\mathbf{p}d\mathbf{q}, $$

with

$$ w^{\epsilon}(\mathbf{q}, \mathbf{p}) = \int_{\mathbb{R}^d} e^{-\frac{1}{\epsilon}\mathbf{p} \cdot (\mathbf{y} - \mathbf{q}) - \frac{1}{2\epsilon}|\mathbf{y} - \mathbf{q}|^2} \psi_0^{\epsilon}(\mathbf{y}) d\mathbf{y}. \qquad (2.2) $$

The Hamiltonian flow is given by

$$ \frac{d\mathbf{Q}}{dt} = \mathbf{P}, \qquad \frac{d\mathbf{P}}{dt} = -\partial_{\mathbf{Q}}V(\mathbf{Q}, t), \qquad (2.3) $$

with initial data

$$ \mathbf{Q}(0, \mathbf{q}, \mathbf{p}) = \mathbf{q}, \qquad \mathbf{P}(0, \mathbf{q}, \mathbf{p}) = \mathbf{p}. \qquad (2.4) $$

The classical action function satisfies

$$ \frac{dS}{dt} = \frac{|\mathbf{P}|^2}{2} - V(\mathbf{Q}, t), \qquad (2.5) $$

and the amplitude $a$ is solution to

$$ \frac{da}{dt}(t, \mathbf{q}, \mathbf{p}) = \frac{1}{2}a(t, \mathbf{q}, \mathbf{p})\text{Tr}(\mathbf{Z}^{-1}(\partial_{\mathbf{z}}\mathbf{P} - i\partial_{\mathbf{z}}\mathbf{Q}\partial_{\mathbf{q}}^2\mathbf{Q})) $$

with

$$ \partial_{\mathbf{z}} = \partial_{\mathbf{q}} - i\partial_{\mathbf{p}}, \qquad \mathbf{Z} = \partial_{\mathbf{z}}(\mathbf{Q} + i\mathbf{P}), \qquad (2.6) $$

where $\partial_{\mathbf{p}}$ (resp. $\partial_{\mathbf{q}}$) denotes the derivative with respect to $\mathbf{p}$ (resp. $\mathbf{q}$).

We now define $\phi_{\text{FGA}}^{\epsilon, \pm}$ using the operator $\mathcal{I}^{\epsilon, \pm}(t; a)$ as follows

$$ \mathcal{I}^{\epsilon, \pm}(t; a)\psi_0^{\epsilon}(\mathbf{x}) := \phi_{\text{FGA}}^{\epsilon, \pm}(t, \mathbf{x}) = \frac{1}{(2\pi\epsilon)^{3d/2}} \int_{\mathbb{R}^d \times \Omega_{\mathbf{q}}^{\pm}} a(t, \mathbf{q}, \mathbf{p})w^{\epsilon}(\mathbf{q}, \mathbf{p})e^{\frac{i}{\epsilon}(S(t, \mathbf{q}, \mathbf{p}) + \mathbf{P}(t, \mathbf{q}, \mathbf{p}) \cdot (\mathbf{x} - \mathbf{Q}(t, \mathbf{q}, \mathbf{p})))} \\ \times e^{-\frac{1}{2\epsilon}|\mathbf{x} - \mathbf{Q}(t, \mathbf{q}, \mathbf{p})|^2} d\mathbf{p}d\mathbf{q}, \qquad (2.7) $$
where $w^\epsilon$ is given by (2.2), and $\Omega_q^\pm$ are two subdomains such that $\mathbb{R}^d = \Omega_q^+ \cup \Omega_q^-$ and $\Omega_q^+ \cap \Omega_q^- = \emptyset$. Naturally, $\phi_{\text{FGA}}^{\epsilon,+}$ and $\phi_{\text{FGA}}^{\epsilon,-}$ can be computed independently at any time. We detail below the overall computational complexity.

**One-domain approach (sequential computing).** From a numerical point of view and following the notations in [43], with indices $(\mathbf{j}, \mathbf{k}, \mathbf{l})$ as elements of index sets denoted by $D_\mathbf{x} \times D_\mathbf{q} \times \tilde{D}_\mathbf{p}^1$, the FGA reads at time $t^n$ and discretization nodes $\mathbf{x}^j$ as

$$ \phi_{\text{FGA}}^{\epsilon}(t^n, \mathbf{x}^j) = \sum_{(\mathbf{k}, \mathbf{l}) \in D_{\mathbf{q}} \times \tilde{D}_{\mathbf{p}}} \frac{a^{\mathbf{k},1}(t^n) r_{\theta}}{(2\pi\epsilon)^{3d/2}} w^{\epsilon}(\mathbf{q}^{\mathbf{k}}, \mathbf{p}^{\mathbf{k},1}) e^{\frac{i}{\epsilon} (S^{\mathbf{k},1}(t^n) + \mathbf{P}^{\mathbf{k},1}(t^n) \cdot (\mathbf{x}^j - \mathbf{Q}^{\mathbf{k},1}(t^n))) - \frac{1}{2\epsilon} |\mathbf{x}^j - \mathbf{Q}^{\mathbf{k},1}(t^n)|^2} |\delta\mathbf{q}| |\delta\mathbf{p}|, \quad (2.8) $$

where $a^{\mathbf{k},1}(t^n)$ (resp. $\mathbf{P}^{\mathbf{k},1}(t^n)$, $\mathbf{Q}^{\mathbf{k},1}(t^n)$) denotes the approximation of $a$ (resp. $\mathbf{P}$, $\mathbf{Q}$) at $(\mathbf{q}^{\mathbf{k}}, \mathbf{p}^{\mathbf{k},1})$ and time $t^n$. Then $(\mathbf{Q}^{\mathbf{k},1}(t^n), \mathbf{P}^{\mathbf{k},1}(t^n))$ and $a^{\mathbf{k},1}(t^n)$ are updated using a fourth order Runge-Kutta (RK4) scheme, and $r_\theta$ is the local truncation function with a radius $\theta$.

The overall computational cost to compute $\phi_{\text{FGA}}^{\epsilon}(t^n, \mathbf{x}^j)$ is estimated as follows. For $(\mathbf{k}, \mathbf{l}) \in D_{\mathbf{q}} \times \tilde{D}_{\mathbf{p}}$, one needs

* $N_{w^\epsilon}^{\mathbf{k},1}$ operations to compute $w^\epsilon(\mathbf{q}^{\mathbf{k}}, \mathbf{p}^{\mathbf{k},1})$;

* $N_{\mathbf{Q}}^{n,\mathbf{k},1}$ RK4 operations to solve $d\mathbf{Q}^{\mathbf{k},1}/dt = \mathbf{P}^{\mathbf{k},1}$ up to time $t^n$;

* $N_{\mathbf{P}}^{n,\mathbf{k},1}$ RK4 operations to solve $d\mathbf{P}^{\mathbf{k},1}/dt = -\partial_{\mathbf{Q}}V(\mathbf{Q}^{\mathbf{k},1})$ up to time $t^n$;

* $N_{\mathbf{a}}^{n,\mathbf{k},1}$ RK4 operations to solve

$$ \frac{da^{\mathbf{k},1}}{dt} = \frac{1}{2} a^{\mathbf{k},1} \text{Tr} \left( (\mathbf{Z}^{-1})^{\mathbf{k},1} (\partial_z \mathbf{P}^{\mathbf{k},1} - i \partial_z \mathbf{Q}^{\mathbf{k},1} \partial_{\mathbf{Q}}^2 V^{\mathbf{k},1}) \right), $$

where $V^{\mathbf{k},1}$ (resp. $(\mathbf{Z}^{-1})^{\mathbf{k},1}$) denotes the value of $V$ (resp. $\mathbf{Z}^{-1}$) at $(\mathbf{q}^{\mathbf{k}}, \mathbf{p}^{\mathbf{k},1})$;

* $cN_{\mathbf{q}}\tilde{N}_{\mathbf{p}}$ operations to compute $\phi_{\text{FGA}}^{\epsilon}(t^n, \mathbf{x}^j)$ by (2.8), where $c$ is a fixed integer, $N_{\mathbf{q}}$ is the cardinal of $D_{\mathbf{q}}$, and $\tilde{N}_{\mathbf{p}}$ is the cardinal of $\tilde{D}_{\mathbf{p}}$.

As a consequence, the total number of operations to compute the FGA up to time $t^n$ and in $\mathbf{x}^j$ is

$$ \text{Op} = \text{Op}(\phi_{\text{FGA}}^{\epsilon}(t^n, \mathbf{x}^j)) \approx \sum_{(\mathbf{k}, \mathbf{l}) \in D_{\mathbf{q}} \times \tilde{D}_{\mathbf{p}}} \left( N_{w^{\epsilon}}^{\mathbf{k},1} + N_{\mathbf{Q}}^{n,\mathbf{k},1} + N_{\mathbf{P}}^{n,\mathbf{k},1} + N_{\mathbf{a}}^{n,\mathbf{k},1} \right) + cN_{\mathbf{q}}\tilde{N}_{\mathbf{p}}. \quad (2.9) $$

**Two-domain approach (parallel computing).** Using the above notations, we directly have in $\Omega_\eta^\pm$

$$ \phi_{\text{FGA}}^{\epsilon, \pm}(t^n, \mathbf{x}^j) = \sum_{(\mathbf{k}, \mathbf{l}) \in D_{\mathbf{q}}^{\pm} \times \tilde{D}_{\mathbf{p}}} \frac{a^{\mathbf{k},1}(t^n) r_{\theta}}{(2\pi\epsilon)^{3d/2}} w^{\epsilon}(\mathbf{q}^{\mathbf{k}}, \mathbf{p}^{\mathbf{k},1}) \times e^{\frac{i}{\epsilon} (S^{\mathbf{k},1}(t^n) + \mathbf{P}^{\mathbf{k},1}(t^n) \cdot (\mathbf{x}^j - \mathbf{Q}^{\mathbf{k},1}(t^n))) - \frac{1}{2\epsilon} |\mathbf{x}^j - \mathbf{Q}^{\mathbf{k},1}(t^n)|^2} |\delta\mathbf{q}| |\delta\mathbf{p}|, \quad (2.10) $$

where for $(\mathbf{k}, \mathbf{l}) \in D_{\mathbf{q}}^{\pm} \times \tilde{D}_{\mathbf{p}}$, we define $N_{w^\epsilon}^{\mathbf{k},1}$, $N_{\mathbf{Q}}^{n,\mathbf{k},1}$, $N_{\mathbf{P}}^{n,\mathbf{k},1}$, and $N_{\mathbf{P}}^{n,\mathbf{k},1}$ as above. As a consequence, the total number of operations to update the FGA up to time $t^n$ and in $\mathbf{x}^j$ is $\text{Op} = \text{Op}^+ + \text{Op}^-$, where

$$ \text{Op}^{\pm} = \text{Op}(\phi_{\text{FGA}}^{\epsilon, \pm}(t^n, \mathbf{x}^j)) \approx \sum_{(\mathbf{k}, \mathbf{l}) \in D_{\mathbf{q}}^{\pm} \times \tilde{D}_{\mathbf{p}}} \left( N_{w^{\epsilon}}^{\mathbf{k},1} + N_{\mathbf{Q}}^{n,\mathbf{k},1} + N_{\mathbf{P}}^{n,\mathbf{k},1} + N_{\mathbf{a}}^{n,\mathbf{k},1} \right) + cN_{\mathbf{q}}^{\pm} \tilde{N}_{\mathbf{p}}, $$

<footnote n="1">¹The notation $\tilde{D}_{\mathbf{p}}$ is used to precise that the integrand in the FGA is localized with respect to $\mathbf{p}$, see [36, 43].</footnote>
where $N_q^\pm$ is the cardinal of $D_q^\pm$.

We observe that $\Omega_\eta^\pm$ can be chosen such that $N_q^+ \approx N_q^-$, so that $\text{Op}^+ \approx \text{Op}^-$ implying that $\text{Op}^+ \approx \text{Op}/2$. We then expect a quasi-ideal speed-up in the case of a parallel implementation of the FGA method. Compared to the standard DDM [12], this method allows for avoiding Schwarz' iterations. FGA may appear as computationally complex since it lives in phase space. However, when the initial data is localized in phase space, the computational complexity for (2.8) can be drastically reduced; see [36, 43]. In addition, we will show that domain decomposition algorithms with ideal speed-up can be easily derived. Numerical examples for this method will be given in Section 6. We show in Section 4 that this DDM can also be useful in the quantum framework to efficiently implement in a parallel framework the use of the FGA-based Semiclassical Schwarz Waveform Relaxation algorithm (4.1).

### 3. Analogies of transparent boundary conditions between FGA and OSWR

In this section, we analyze the analogy of TBCs between the formalisms of FGA [43] and OSWR [12]. This is a preparation for the development of Semiclassical Schwarz Waveform Relaxation (SSWR) methods in the later two sections. Recall that OSWR algorithms are derived from TBCs or high order ABCs used to avoid spurious reflections at the boundaries [12]. We shall show that SSWR are actually derived from OSWR, using similar Hamiltonian flow and outgoing wave concepts in the semiclassical limit.

We denote by $\Omega$ an open bounded domain in $\mathbb{R}^d$. We show analogies between TBCs on $\partial\Omega$ in the FGA formalism in the semiclassical regime [25, 39] and TBCs derived using microlocal analysis beyond the semiclassical regime [1, 7]. We introduce a set $\Omega^\varepsilon$ containing $\Omega$, defined as $\Omega^\varepsilon = \cup_{x\in\Omega} B(x, \sqrt{\varepsilon})$, where $B(x, \sqrt{\varepsilon})$ is a ball of center $x$ and radius $\sqrt{\varepsilon}$. We recall that the FGA approximating the solution to (1.1) reads as

$$ \phi_{\text{FGA}}^{\varepsilon}(t, \mathbf{x}) = \frac{1}{(2\pi\varepsilon)^{3d/2}} \int_{\mathbb{R}^{3d}} a(t, \mathbf{q}, \mathbf{p}) e^{\frac{i}{\varepsilon} \Phi(t, \mathbf{x}, \mathbf{y}, \mathbf{q}, \mathbf{p})} \phi_0(\mathbf{y}) d\mathbf{y} d\mathbf{p} d\mathbf{q}, $$

with the phase function $\Phi$ given by

$$ \begin{aligned} \Phi(t, \mathbf{x}, \mathbf{y}, \mathbf{q}, \mathbf{p}) ={}& S(t, \mathbf{q}, \mathbf{p}) + \frac{i}{2} |\mathbf{x} - \mathbf{Q}(t, \mathbf{q}, \mathbf{p})|^2 + \mathbf{P}(t, \mathbf{q}, \mathbf{p}) \cdot (\mathbf{x} - \mathbf{Q}(t, \mathbf{q}, \mathbf{p})) \\ & + \frac{i}{2} |\mathbf{y} - \mathbf{q}|^2 - \mathbf{p} \cdot (\mathbf{y} - \mathbf{q}), \end{aligned} \quad (3.1) $$

where $\mathbf{P}, \mathbf{Q}, \mathbf{S}$ are respectively defined in (2.3), (2.4), (2.5). In [43], outgoing wave conditions at time $t$, for $B(\mathbf{Q}_0(t, \mathbf{q}, \mathbf{p}), \sqrt{\varepsilon}) \in \Omega_i^\varepsilon - \Omega$, are determined according to the condition

$$ \frac{d\mathbf{Q}_0}{dt} \cdot \mathbf{n}_0 > 0, \quad (3.2) $$

where $\mathbf{n}_0$ is the outward normal vector of $\partial\Omega$ at $\mathbf{Q}(t, \mathbf{q}, \mathbf{p})$. In that case, the contribution of $\mathbf{Q}_0$ is removed from (3.1).

Let us now define the Schrödinger operator by

$$ P_{\varepsilon} = -i\varepsilon\partial_t - \frac{\varepsilon^2}{2}\Delta_{\mathbf{x}} + V(\mathbf{x}, t). \quad (3.3) $$

Then, its classical symbol is: $p_\varepsilon(\mathbf{x}, t, \xi, \omega) = \frac{\varepsilon^2}{2} |\xi|^2 + \varepsilon\omega + V(\mathbf{x}, t)$. We denote by $T^*([0, T] \times \partial\Omega)$ the cotangent bundle such that the rays propagate in the hyperbolic zone [27, 41]

$$ \mathcal{H} = \{(\mathbf{x}_0, t_0, \xi_0, \omega_0) \in T^*([0, T] \times \partial\Omega) : \varepsilon\omega_0 + \varepsilon^2 |\xi_0|^2/2 + V < 0\}. \quad (3.4) $$
The outgoing bicharacteristic strips (from which the TBCs (2.3) are derived) are then selected from the following Hamiltonian equations ([26]), for $s \in \mathbb{R}$,

$$ egin{cases} \frac{d\mathbf{x}}{ds}(s) = \partial_{\xi} p_{\varepsilon}(\mathbf{x}(s), t(s), \xi(s), \omega(s)) = \varepsilon^2 \xi(s), \\ \frac{dt}{ds}(s) = \partial_{\omega} p_{\varepsilon}(\mathbf{x}(s), t(s), \xi(s), \omega(s)) = \varepsilon, \\ \frac{d\xi}{ds}(s) = -\partial_{\mathbf{x}} p_{\varepsilon}(\mathbf{x}(s), t(s), \xi(s), \omega(s)) = -\partial_{\mathbf{x}} V(\mathbf{x}, t(s)), \\ \frac{d\omega}{ds}(s) = -\partial_{t} p_{\varepsilon}(\mathbf{x}(s), t(s), \xi(s), \omega(s)) = 0, \end{cases} \quad (3.5) $$

based on the condition at $\mathbf{x}_0 \in \partial\Omega$, with outward normal vector $\mathbf{n}_0$ at $\partial\Omega$,

$$ \frac{d\mathbf{x}}{ds}(s) \cdot \mathbf{n}_0 > 0. \quad (3.6) $$

Note that $\omega(s) = \omega_0$ and $t(s) = \varepsilon s + t_0$, and an analogy can now be deduced between (3.2) and (3.6). Indeed, from (2.3)

$$ \mathbf{P} = \frac{d\mathbf{Q}}{dt} = \frac{d\mathbf{Q}}{ds}\frac{ds}{dt} = \frac{1}{\varepsilon}\frac{d\mathbf{Q}}{ds}, \quad -\partial_{\mathbf{Q}}V(\mathbf{Q}, t) = \frac{d\mathbf{P}}{dt} = \frac{d\mathbf{P}}{ds}\frac{ds}{dt} = \frac{1}{\varepsilon}\frac{d\mathbf{P}}{ds}, $$

one has

$$ \frac{d\mathbf{Q}}{ds} = \varepsilon\mathbf{P}, \quad \frac{d\mathbf{P}}{ds} = -\varepsilon\partial_{\mathbf{Q}}V(\mathbf{Q}, t), $$

which exhibits an analogy between $\mathbf{x}(s)$ and $\mathbf{Q}(s)$, and between $\xi(s)$ and $\mathbf{P}(s)/\varepsilon$.

As in [4], the Schrödinger operator can be factorized, in the hyperbolic region of the cotangent bundle $T^*([0, T] \times \partial\Omega)$ given by (3.4), as the composition of an outgoing and an incoming wave operators at boundaries, and specifically in one dimension,

$$ i\varepsilon\partial_t + \frac{\varepsilon^2}{2}\partial_x^2 - V = \left(\frac{\varepsilon}{\sqrt{2}}\partial_x + \sqrt{V - i\varepsilon\partial_t}\right)\left(\frac{\varepsilon}{\sqrt{2}}\partial_x - \sqrt{V - i\varepsilon\partial_t}\right). \quad (3.7) $$

Then, the OSWR algorithm [12, 24] can be derived from TBCs based on DtN operators by using the above factorization

$$ \left\{ \begin{array}{l} i\varepsilon\partial_t \phi^{\varepsilon, \pm, (k)} \\ \left( \frac{\varepsilon}{\sqrt{2}}\partial_x \pm \sqrt{V - i\varepsilon\partial_t} \right) \phi^{\varepsilon, \pm, (k)} \end{array} \right. = \begin{array}{l} H_{\varepsilon}\phi^{\varepsilon, \pm, (k)} \quad \text{in } (0, T) \times \Omega_{\eta}^{\pm}, \\ \left( \frac{\varepsilon}{\sqrt{2}}\partial_x \pm \sqrt{V - i\varepsilon\partial_t} \right) \phi^{\varepsilon, \mp, (k-1)} \quad \text{on } (0, T) \times \Gamma_{\eta}^{\pm}. \end{array} \quad (3.8) $$

Therefore, according to the above discussions, the key of developing semiclassical Schwarz Waveform Relaxation method lies in mimicking OSWR algorithms using FGA-based conditions rather than DtN operators.

## 4. Herman-Kluk Schwarz Waveform Relaxation algorithm in the middle-frequency regime

In this and the next sections, we systematically discuss the development of semiclassical Schwarz Waveform Relaxation (SWR) methods for solving (1.1) in the middle frequency regime. Based on the analysis given in Section 3, a natural way is to approximate the transparent boundary conditions given in (3.8) by frozen Gaussian approximation, which yields the following algorithm.
### 4.1. A first FGA-based domain decomposition algorithm in the middle-frequency regime

Using the notations introduced in Section 2 and for $x_\eta \in \Gamma_\eta^\pm$, the FGA reads

$$ \mathcal{I}^{\varepsilon}(t; a)\phi_{0}^{\varepsilon}(x_{\eta}^{\pm}) = \frac{1}{(2\pi\varepsilon)^{3d/2}} \int_{\mathbb{R}^{3d}} a(t, \mathbf{q}, \mathbf{p})e^{\frac{1}{\varepsilon}\Phi(t, x_{\eta}^{\pm}, \mathbf{y}, \mathbf{q}, \mathbf{p})} \phi_{0}^{\varepsilon}(\mathbf{y})d\mathbf{y}d\mathbf{p}d\mathbf{q}, $$

with $\Phi$ given in (3.1). For simplicity of notations, we shall omit hereafter the index $\varepsilon$ in $\phi$ without causing any confusion. Then the FGA-based SWR algorithm can be designed as follows

$$ \left\{ \begin{array}{ll} \dot{i}\varepsilon\partial_t\phi^{\pm,(k)}(t, \cdot) & = H_{\varepsilon}\phi^{\pm,(k)}(t, \cdot) \text{ in } \{t > t^n\} \times \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(t, \cdot) & = \mathcal{I}^{\varepsilon}(t; a)\phi^{(k-1)}(t^n, \cdot) \text{ on } \{t > t^n\} \times \Gamma_{\eta}^{\pm}, \end{array} \right. \qquad (4.1) $$

where we assumed that, at Schwarz iteration $k$ and time $t^n$, the solution $\phi^{\pm,(k)}(t^n, \cdot)$ is known in $\Omega_\eta^\pm$, as well as $S(t^n, \mathbf{q}, \mathbf{p})$, $\mathbf{Q}(t^n, \mathbf{q}, \mathbf{p})$, $\mathbf{P}(t^n, \mathbf{q}, \mathbf{p})$, $a(t^n, \mathbf{q}, \mathbf{p})$, with $(\mathbf{q}, \mathbf{p}) \in \Omega_\mathbf{q} \times \Omega_\mathbf{p}$. In (4.1), we denoted by $\phi^{(k-1)}$ the reconstructed global solution in $\mathbb{R}^d$ at Schwarz iteration $(k-1)$, which is typically constructed as follows

$$ \phi^{(k-1)}(t, \mathbf{y}) = \begin{cases} \frac{\phi^{\pm,(k-1)}(t, \mathbf{y})}{\phi^{+, (k-1)} + \phi^{-, (k-1)}} & \text{for } (t, \mathbf{y}) \in (t_n, t) \times \Omega^{\pm} \setminus (\Omega^+ \cap \Omega^-), \\ \frac{\phi^{+, (k-1)} + \phi^{-, (k-1)}}{2} & \text{for } (t, \mathbf{y}) \in (t_n, t) \times \Omega^+ \cap \Omega^- \end{cases} \qquad (4.2) $$

Then, for $x_\eta^\pm \in \Gamma^\pm$ and $t > t^n$, we have

$$ \mathcal{I}^{\varepsilon}(t; a)\phi^{(k-1)}(t^n, x_{\eta}^{\pm}) = \frac{1}{(2\pi\varepsilon)^{3d/2}} \int_{\mathbb{R}^{2d}} a(t, \mathbf{q}, \mathbf{p})w^{\varepsilon,(k-1)}(t^n, \mathbf{q}, \mathbf{p}) \times e^{\frac{1}{\varepsilon}\left(S(t, \mathbf{q}, \mathbf{p}) + \mathbf{P} \cdot (x_{\eta}^{\pm} - \mathbf{Q}(t, \mathbf{q}, \mathbf{p}))\right)} \times e^{-\frac{1}{2\varepsilon}|x_{\eta}^{\pm} - \mathbf{Q}(t, \mathbf{q}, \mathbf{p})|^2} d\mathbf{p}d\mathbf{q}, \qquad (4.3) $$

where

$$ w^{\varepsilon,(k-1)}(t^n, \mathbf{q}, \mathbf{p}) := \int_{\Omega_{\eta}^{\pm}} \exp \left( -\frac{\dot{i}}{\varepsilon} \mathbf{p} \cdot (\mathbf{y} - \mathbf{q}) - \frac{1}{2\varepsilon} |\mathbf{y} - \mathbf{q}|^2 \right) \phi^{(k-1)}(t^n, \mathbf{y}) d\mathbf{y}. $$

We discretize the transmission condition (4.3) as follows. Let us denote $(\mathbf{k}, \mathbf{l}) \in D_\mathbf{q} \times \tilde{D}_\mathbf{p}$ as the discretization indices in $(\mathbf{q}, \mathbf{p}) \in \Omega_\mathbf{q} \times \Omega_\mathbf{p}$ for some $t^{n+1} > t^n$, and $\mathbf{j} \in D_\mathbf{y}$ as spatial indices for $\mathbb{R}^d$. We then have

$$ \mathcal{I}^{\varepsilon}(t^{n+1}; a)\phi^{(k-1)}(t^n, x_{\eta}^{\pm}) \approx \frac{1}{(2\pi\varepsilon)^{3d/2}} \sum_{\mathbf{j} \in D_{\mathbf{x}}} \phi_{\mathbf{j}}^{(k-1)}(t^n) \times \sum_{(\mathbf{k}, \mathbf{l}) \in D_{\mathbf{q}} \times \tilde{D}_{\mathbf{p}}} \alpha_{\mathbf{j}}(\mathbf{q}^{\mathbf{k}}, \mathbf{p}^{\mathbf{k},\mathbf{l}})e^{\frac{1}{\varepsilon}\left(S^{n+1,\mathbf{k},\mathbf{l}} + \mathbf{P}^{n+1,\mathbf{k},\mathbf{l}} \cdot (x_{\eta}^{\pm} - \mathbf{Q}^{n+1,\mathbf{k},\mathbf{l}})\right)} \times e^{-\frac{1}{2\varepsilon}|x_{\eta}^{\pm} - \mathbf{Q}^{n+1,\mathbf{k},\mathbf{l}}|^2} |\delta\mathbf{q}| |\delta\mathbf{p}|, $$

where

$$ \alpha_{\mathbf{j}}(\mathbf{q}^{\mathbf{k}}, \mathbf{p}^{\mathbf{k},\mathbf{l}}) = |\delta\mathbf{y}|e^{-\frac{1}{\varepsilon}\mathbf{p}^{\mathbf{k},\mathbf{l}} \cdot (\mathbf{y}_{\mathbf{j}} - \mathbf{q}^{\mathbf{k}}) - \frac{1}{2\varepsilon}|\mathbf{y}_{\mathbf{j}} - \mathbf{q}^{\mathbf{k}}|^2}. $$

This integral can easily be performed in parallel (in $\mathbf{q}$) as proposed in Section 2. Formally, the transmission conditions read

$$ \phi^{\pm,(k)}(t^{n+1}, x_{\eta}^{\pm}) \approx \sum_{\mathbf{j}} \beta_{\mathbf{j}}^{\pm}(t^n) \phi_{\mathbf{j}}^{(k-1)}(t^n), $$
with $\beta_j^{\pm}$ defined as

$$ \beta_j^{\pm}(t^n) = \frac{1}{(2\pi\epsilon)^{3d/2}} \sum_{(k,l) \in D_q \times \tilde{D}_p} \alpha_j(\mathbf{q}^k, \mathbf{p}^{k,l}) e^{\frac{i}{\epsilon} (S^{n+1,k,l} + \mathbf{p}^{n+1,k,l} \cdot (\mathbf{x}_{\eta}^{\pm} - \mathbf{Q}^{n+1,k,l}))} - \frac{1}{2\epsilon} |\mathbf{x}_{\eta}^{\pm} - \mathbf{Q}^{n+1,k,l}|^2 |\delta\mathbf{q}| |\delta\mathbf{p}|. $$

The above approach comes straightforwardly from the analogies given in Section 3 and is easy to be implemented. However, it suffers from an issue of non-convergence as analyzed below.

**Proposition 4.1.** We denote the Schrödinger operator by $P_{\epsilon}$ as in (3.3) and we assume that $V$ is a sub-quadratic potential in $C^{\infty}(\mathbb{R}^d, \mathbb{R})$. Then, for any $\phi_0 \in L^2(\mathbb{R}^d, \mathbb{C})$, the Schwarz algorithm, for $k \ge 1$,

$$ \left\{ \begin{array}{ll} P_{\epsilon}\phi^{\pm,(k)}(t, \cdot) & = 0, & \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(0, \cdot) & = \phi_0(\cdot), & \text{on } \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(t, \cdot) & = \mathcal{I}^{\epsilon}(t; a)\phi_0, & \text{on } (0, T) \times \Gamma_{\eta}^{\pm}, \end{array} \right. \quad (4.4) $$

is not convergent.

**Proof.** Recall that $\mathcal{I}^{\epsilon}(t; a)\phi_0$ is the FGA solution. From [39], one has at any time $T > 0$,

$$ \sup_{t \in (0, T)} \left\| \Psi_e(t, \cdot) - \mathcal{I}^{\epsilon}(t; a) \cdot \phi_0 \right\|_{L^2(\mathbb{R}^d)} \le c(T) \epsilon \|\phi_0\|_{L^2(\mathbb{R}^d)}, $$

where we denote by $\Psi_e$ the exact solution. At any time $t \in (0, T)$ and for any $\mathbf{x}_{\eta}^{\pm} \in \Gamma_{\eta}^{\pm}$, we have

$$ \Psi_e(t, \mathbf{x}_{\eta}^{\pm}) - \mathcal{I}^{\epsilon}(t; a)\phi_0(\mathbf{x}_{\eta}^{\pm}) = \epsilon c_{\eta}(t, \mathbf{x}_{\eta}^{\pm}), \quad (4.5) $$

where $c_{\eta}(t, \cdot) \in L^2(\mathbb{R}^d)$.

We define $e^{\pm,(k)} := \Psi_e|_{\Omega_{\eta}^{\pm}} - \phi^{\pm,(k)}$, and then

$$ \left\{ \begin{array}{ll} P_{\epsilon}e^{\pm,(k)}(t, \cdot) & = 0, & \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ e^{\pm,(k)}(0, \cdot) & = 0, & \text{on } \Omega_{\eta}^{\pm}, \\ e^{\pm,(k)}(t, \cdot) & = \epsilon c_{\eta}(t, \cdot), & \text{on } (0, T) \times \Gamma_{\eta}^{\pm}. \end{array} \right. $$

The error is non-null and $k$-independent and the algorithm is then non-convergent. $\square$

We remark that, although this method is not convergent, the error is expected to be negligible for small values of $\epsilon$ and $\eta$, which will be numerically tested in Section 6.

## 4.2. Herman-Kluk Schwarz Waveform Relaxation methods

In order to derive a convergent Schwarz method, we need to introduce a $k$-dependence in the transmission conditions. In addition as, on coarse grids, FGA computes (less oscillatory) amplitudes more accurate than phases of high oscillatory wavefunctions, we propose to use the modulus of FGA only and refer the method to as Herman-Kluk Schwarz Waveform Relaxation methods (HKSWR).

We first introduce a function $\mathcal{F}_{\gamma}: \mathbb{C} \times \mathbb{C} \to \mathbb{R}$, $\gamma \in [0, 1]$, $\forall (z_1, z_2) \in \mathbb{C}^2$, defined by

$$ \mathcal{F}_{\gamma}(z_1, z_2) = \begin{cases} z_2 \left| \frac{z_1}{z_2} \right|^{\gamma}, & \text{if } z_2 \neq 0, \\ 0, & \text{if } z_2 = 0. \end{cases} \quad (4.6) $$

Note that for any $\gamma = 1$, $\mathcal{F}_1(z_1, z_2) = |z_1| \exp(\dot{\mathbf{i}} \arg(z_2))$.

We consider the following iterative scheme, for $k \ge 1$ and any initial data $\phi_0 \in L^2(\mathbb{R}^d, \mathbb{C})$,

$$ \left\{ \begin{array}{ll} P_{\epsilon}\phi^{\pm,(k)} & = 0, & \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(0, \cdot) & = \phi_0(\cdot), & \text{on } \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(t, \cdot) & = \mathcal{F}_{\gamma_{\epsilon}(k)}(\mathcal{I}^{\epsilon}(t; a)\phi_0(\cdot), \phi^{(k-1)}(t, \cdot)), & \text{on } (0, T) \times \Gamma_{\eta}^{\pm}, \end{array} \right. \quad (4.7) $$
with $P_\varepsilon$ defined in (3.3) and $\phi^{(k-1)}$ in (4.2). Here $\gamma_\varepsilon$ is a function of $k$: $\mathbb{N}^* \mapsto [0, 1]$, which shall be specified later in Section 4.3 for fixed $\varepsilon > 0$. Remark that choosing the initial guess as $\phi^{\pm,(0)}(t, \cdot) = \mathcal{I}^\varepsilon(t; a)\phi_0$ and $\gamma = 1$ will degenerate the algorithm back to (4.4), and thus (4.7) is not convergent if $\gamma_\varepsilon$ is identically equal to one, according to Proposition 4.1. According to [35], the smaller $\varepsilon$ the more accurate FGA. Thus we expect that the contribution from FGA will lead to a small error after a few Schwarz iterations. To ensure the convergence of the algorithm, when $k$ gets large, one needs to switch to the CSWR algorithm, which reads as follows

$$ \left\{ \begin{array}{lll} P_{\varepsilon}\phi^{\pm,(k)}(t, \cdot) & = & 0, \quad \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(0, \cdot) & = & \phi_0(\cdot), \quad \text{on } \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(t, \cdot) & = & \phi^{\mp,(k-1)}(t, \cdot), \quad \text{on } (0, T) \times \Gamma_{\eta}^{\pm}. \end{array} \right. \qquad (4.8) $$

The error equation reads for $k \ge 2$

$$ \left\{ \begin{array}{lll} P_{\varepsilon}e^{\pm,(k)}(t, \cdot) & = & 0, \quad \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ e^{\pm,(k)}(0, \cdot) & = & 0, \quad \text{on } \Omega_{\eta}^{\pm}, \\ e^{\pm,(k)}(t, \cdot) & = & e^{\pm,(k-2)}(t, \cdot), \quad \text{on } (0, T) \times \Gamma_{\eta}^{\pm}, \end{array} \right. $$

and the convergence of CSWR algorithm (4.8) is assured by the following theorem, with proofs in [22]; see also [11] for the non-constant potential case.

**Theorem 4.1.** Assume that $V$ is a constant potential. Starting from a null initial guess in $\Omega_\eta^\pm$, the CSWR algorithm (4.8) is convergent and the following inequality holds

$$ \|(e^{+, (2k+1)}, e^{-, (2k+1)})\|_{H^{3,3/2}(\Omega_{\eta}^{+} \times (0,T)) \times H^{3,3/2}(\Omega_{\eta}^{-} \times (0,T))} \leq D C_{\varepsilon, \eta}^{k} \|(e^{+, (0)}, e^{-, (0)})\|_{(H^{3/4}(0,T))^{2}}, \quad (4.9) $$

where $D$ is a real-valued positive constant and coefficient $C_{\varepsilon,\eta}$ is a contraction factor, which is defined by $C_{\varepsilon,\eta} = \sup_{\tau \in \mathcal{H}_\tau} |L_{\varepsilon,\eta}(\tau)|$, for $|\tau| \gg 1$, where $\mathcal{H}_\tau = \{\tau : V + \varepsilon\tau + \varepsilon^2\xi^2 < 0\}$ and

$$ L_{\varepsilon, \eta}(\tau) \approx \exp \left( -\eta \sqrt{-2\tau/\varepsilon} - V\eta \frac{1}{\sqrt{-2\varepsilon^3\tau}} \right). \qquad (4.10) $$

Therefore, it is natural to choose $\gamma_\varepsilon(k)$ close to 1 for small $k$, since Proposition 4.1 shows that, when $\varepsilon$ is small, the reconstructed solution will be close (but not equal) to the exact one after the first few Schwarz iterations. For larger $k$, Theorem 4.1 suggests to take $\gamma_\varepsilon$ close to 0 to ensure the convergence.

**Remark 4.1.** The transmission conditions of HKSWR can also be chosen of the form

$$ \phi^{\pm,(k)}(t, \cdot) = [\mathcal{I}^{\varepsilon}(t; a)\phi_0(\cdot)]^{\gamma_{\varepsilon}(k)} [\phi^{\mp,(k-1)}(t, \cdot)]^{1-\gamma_{\varepsilon}(k)}. $$

However, in this case, one needs a fine grid for a discrete accurate representation of the FGA solution due to the highly oscillatory phase function.

In the following, we will consider a slightly modified version of (4.7)

$$ \left\{ \begin{array}{lll} P_{\varepsilon}\phi^{\pm,(k)}(t, \cdot) & = & 0, \quad \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(0, \cdot) & = & \phi_0(\cdot), \quad \text{on } \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(t, \cdot) & = & \mathcal{F}_{\gamma_{\varepsilon}(k)}(\mathcal{I}^{\tilde{\varepsilon}}(t; a) \cdot \phi_0(\cdot), \phi^{\mp,(k-1)}(t, \cdot)), \quad \text{on } (0, T) \times \Gamma_{\eta}^{\pm}, \end{array} \right. \qquad (4.11) $$

where $\tilde{\varepsilon}$ is a small parameter, which in principle satisfies $\tilde{\varepsilon} \le \varepsilon$. Let us notice that $\tilde{\varepsilon}$ is used as a free parameter in the transmission conditions. In practice, $\tilde{\varepsilon}$ does not need to be equal to $\varepsilon$ since the convergence is finally ensured by the CSWR method. However, the more accurate the FGA, the better (say for $\tilde{\varepsilon} \approx \varepsilon \ll 1$). The ideal choice for $\tilde{\varepsilon}$, in combination with the function $\gamma_\varepsilon$, is then the one that lead to the fastest convergence of the DD algorithm.
### 4.3. Criteria for the $\gamma_{\varepsilon}$-function

This subsection is devoted to the construction of the function $\gamma_{\varepsilon}$ allowing an efficient coupling between the FGA and CSWR transmission conditions, as discussed in Theorem 4.2. First, let us recall that the algorithm based on pure FGA conditions (4.4) is not convergent due to the existence of asymptotic errors of FGA. Theorem 4.1 implies that one can estimate the number of CSWR iterations $k_{\varepsilon,\eta}$ to reach a given error threshold $e_T$

$$k_{\varepsilon,\eta} \approx \frac{\log(e_T/E_0)}{2 \log C_{\varepsilon,\eta}} = \frac{\varepsilon\sqrt{2}}{\eta(V+2\varepsilon)} \log\left(\frac{E_0}{e_T}\right).$$

Then we compare the rate of convergence of HKSWR for different choices of $\gamma_{\varepsilon}$ to validate the above analysis and help to carefully select $\gamma_{\varepsilon}$: i) $\gamma_{\varepsilon}(k) = 0$, for all $k \ge 1$ (CSWR), ii) $\gamma_{\varepsilon}(k) = (-\alpha_{\varepsilon}(k-1)^2)$, iii) $\gamma_{\varepsilon}(k) = 1$, for all $k \ge 1$ (pure FGA), iv) $\gamma_{\varepsilon}(k) = 1$, $\forall k \le k_0$ & $\gamma_{\varepsilon}(k) = 0$, $\forall k > k_0$, with plots of $\gamma_{\varepsilon}$ in Figure 4.1.

<figure><img src="image_3.png" /><figcaption>Figure 4.1: Various functions $\gamma_{\varepsilon}$.</figcaption></figure>

The numerical framework is the same as in Section 6.2, with in addition

* Test 1: $\varepsilon = 2^{-3}$, $\tilde{\varepsilon} = 2^{-6}$, $\alpha_{\varepsilon} = 5 \times 10^{-2}$ and $k_0 = 1$ and 5. The other numerical data are the same as Subsection 6.2.1.

* Test 2: $\varepsilon = 2^{-5}$, $\tilde{\varepsilon} = 2^{-10}$, $\alpha_{\varepsilon} = 10^{-2}$, and $k_0 = 1$ and 10. The other numerical data are the same as Subsection 6.2.2.

The numerical results in Figure 4.2 suggest that smaller $\alpha_{\varepsilon}$ should be chosen at higher frequency (corresponding to smaller $\varepsilon$). This also illustrates the fact that the higher the frequency, the more precise the FGA, and the closer $\gamma_{\varepsilon}$ should be taken to 1. Note also that, as expected, the full FGA approach ($\gamma_{\varepsilon}(k) = 1$ for all $k$) is shown to be numerically non-convergent.

To summarize, one shall construct $\gamma_{\varepsilon}$ according to the following rules: i) $\gamma_{\varepsilon}$ is decreasing in $k$, ii) $\gamma_{\varepsilon}(1) = 1$, iii) and for fixed $\varepsilon$, $\lim_{k \to \infty} \gamma_{\varepsilon}(k) = 0$ and $\sum_{k \in \mathbb{N}^*} \gamma_{\varepsilon}(k) < +\infty$ or such that for some $k_0$, $\gamma_{\varepsilon}(k) = 0$, for all $k \ge k_0$, ensuring the convergence of the HKSWR method, and a possible choice is as follows

$$\gamma_{\varepsilon}(k) = \exp(-\alpha_{\varepsilon}(k-1)^2),$$

where $\alpha_{\varepsilon} \in [0, 1)$ is an increasing function of $\varepsilon$ and goes to 0 when $\varepsilon$ tends toward 0.
<figure><img src="image_4.png" /><figcaption>Figure 4.2: Comparison of the rate of convergence with i) $\varepsilon = 2^{-3}$, $\tilde{\varepsilon} = 2^{-6}$, $\alpha_{\varepsilon} = 5 \times 10^{-2}$, and $k_0 = 1$ and $k_0 = 5$ and ii) $\varepsilon = 2^{-6}$, $\tilde{\varepsilon} = 2^{-10}$, $\alpha_{\varepsilon} = 10^{-2}$, and $k_0 = 1$ and $k_0 = 10$.</figcaption></figure>

### 4.4. Analysis of convergence

In the following theorem, we analyze the convergence of HKSWR in three generic cases.

**Theorem 4.2.** Consider a function $\gamma_{\varepsilon}: \mathbb{N}^* \to [0, 1]$, where is $\varepsilon$ the rescaled Planck constant in (1.1). Then the HKSWR algorithm (4.7) is convergent in the following three generic and distinct cases.

1. Assume that $\gamma_{\varepsilon}$ is chosen such that for some $k_0 \in \mathbb{N}^*$, $\gamma_{\varepsilon}(k) = 0$ for all $k \ge k_0$, then (4.7) is asymptotically convergent as CSWR (4.8).

2. Assume that $\gamma_{\varepsilon}$ is chosen such that for some $k_0 \in \mathbb{N} - \{0, 1\}$, $\gamma_{\varepsilon}(k) = 1$ for all $k \le k_0$, and $\gamma_{\varepsilon}(k) = 0$ for all $k > k_0$, then (4.7) is asymptotically convergent as (4.8), and for $\varepsilon$ small enough, any given tolerance is reached in less iterations than (4.8).

3. Assume that i) $\gamma_{\varepsilon}$ is a decreasing function which is convergent to zero, such that ii) $\sum_{k=1}^{\infty} \gamma_{\varepsilon}(k)$ exists, and that iii) at $(0, T) \times \Gamma_{\eta}^{\pm}$ for all $k \in \mathbb{N}^*$,

$$|e^{\pm,(2k)}| \le |\mathcal{I}^{\varepsilon}(t; a)\phi_0 - \Psi_{\varepsilon}(t, \cdot)|,$$

where $e^{\pm,(2k)}$ is defined in (4.14), then (4.7) is asymptotically convergent as CSWR, and any given tolerance is reached in less iterations than in Case 2.

**Proof.** We successively study the three distinct cases.

*Case 1.* It is easy to see that (4.7) degenerates as above into CSWR for all $k \ge k_0$, and which is known to be convergent [24] independently of the initial guess, with a contraction factor given by (4.10).

*Case 2.* For any $k$ such that $2 \le k \le k_0$, one has

$$ \left\{ \begin{array}{ll} P_{\varepsilon}\phi^{\pm,(k)}(t, \cdot) & = 0, & \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(0, \cdot) & = \phi_0, & \text{on } \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(t, \cdot) & = \mathcal{F}_1(\mathcal{I}^{\varepsilon}(t; a)\phi_0(\cdot), \phi^{\pm,(k-2)}(t, \cdot)), & \text{on } (0, T) \times \Gamma_{\eta}^{\pm}, \end{array} \right. $$

with $\phi^{\pm,(0)}$ given at $\Gamma_{\eta}^{\pm}$ by $\mathcal{I}^{\varepsilon}(t; a)\phi_0(\pm \eta/2)$, which by induction yields

$$ \left\{ \begin{array}{ll} P_{\varepsilon}e^{\pm,(k_0)}(t, \cdot) & = 0, & \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ e^{\pm,(k_0)}(0, \cdot) & = 0, & \text{on } \Omega_{\eta}^{\pm}, \\ e^{\pm,(k_0)}(t, \cdot) & = c_{\eta}(t, \cdot)\varepsilon, & \text{on } (0, T) \times \Gamma_{\eta}^{\pm}, \end{array} \right. $$
where $c_η$ is defined in (4.5). The first $k_0$ iterations allow for an error estimate as small as $\epsilon$. For $k > k_0$, since $\gamma_\epsilon(k) = 0$, HKSWR degenerates into CSWR, and thus Theorem 4.1 gives, for $k > k_0$, that there exists $D > 0$ such that

$$ \|(e^{+, (2k+1)}, e^{-, (2k+1)})\|_{H^{3,3/2}(\Omega_{\eta}^{+} \times (0,T)) \times H^{3,3/2}(\Omega_{\eta}^{-} \times (0,T))} \le \epsilon D C_{\epsilon, \eta}^{k-k_0} \|(c_{\eta}(\cdot, \eta/2), c_{\eta}(\cdot, -\eta/2))\|_{(H^{3/4}(0,T))^2}, $$

which implies that any given tolerance $\delta$ will be reached in less iterations than CSWR when $\epsilon$ is small. The initial use of FGA allows then for a reduction by a $O(\epsilon)$ of the convergence rate, compared to CSWR. Case 3. For $k \ge 1$, we have at $(0, T) \times \Gamma_\eta^\pm$

$$ \phi^{\pm, (k)} = \phi^{\pm, (k-1)} \left| \frac{\mathcal{I}^{\epsilon}(t; a) \phi_0}{\phi^{\pm, (k-1)}} \right|^{\gamma_{\epsilon}(k-1)}. $$

We then have at $(0, T) \times \Gamma_\eta^\pm$ and for $k \ge 2$

$$ \phi^{\pm, (k)} = \phi^{\pm, (k-2)} \left| \frac{\mathcal{I}^{\epsilon}(t; a) \phi_0}{\phi^{\pm, (k-2)}} \right|^{-\gamma_{\epsilon}(k)\gamma_{\epsilon}(k-1)} \left| \mathcal{I}^{\epsilon}(t; a) \phi_0 \right|^{\gamma_{\epsilon}(k)+\gamma_{\epsilon}(k-1)}. \quad (4.12) $$

It turns out that the right hand side in (4.12) can be rewritten in terms of $\gamma_\epsilon$ and $\mathcal{I}^\epsilon(t; a)\phi_0$ only. Indeed, starting from $\phi^{\pm, (0)}(t, \cdot)$ computed from $\mathcal{I}^\epsilon(t; a)\phi_0$, one has

$$ \phi^{\pm, (2)} = \phi^{\pm, (0)} |\mathcal{I}^{\epsilon}(t; a) \phi_0|^{\gamma_{\epsilon}(1)+\gamma_{\epsilon}(2)}. $$

Then since at $(0, T) \times \{\pm \eta/2\}$, $|\phi^{\pm, (0)}| = |\mathcal{I}^{\epsilon}(t; a)\phi_0|$, we still have at $(0, T) \times \Gamma_\eta^\pm$

$$ \phi^{\pm, (4)} = \phi^{\pm, (0)} |\mathcal{I}^{\epsilon}(t; a) \phi_0|^{\sum_{i=1}^{4} \gamma_{\epsilon}(i) + (\gamma_{\epsilon}(1) + \gamma_{\epsilon}(2)) \gamma_{\epsilon}(3) \gamma_{\epsilon}(4)} = \phi^{\pm, (0)} |\mathcal{I}^{\epsilon}(t; a) \phi_0|^{\sum_{i=1}^{4} \gamma_{\epsilon}(i) + \sum_{l=3}^{4} \gamma_{\epsilon}(l) \sum_{i=1}^{2} \gamma_{\epsilon}(i)}, $$

$$ \phi^{\pm, (6)} = \phi^{\pm, (0)} |\mathcal{I}^{\epsilon}(t; a) \phi_0|^{\sum_{i=1}^{6} \gamma_{\epsilon}(i) + \sum_{l=5}^{6} \gamma_{\epsilon}(l) \sum_{i=1}^{4} \gamma_{\epsilon}(i) + (\sum_{l=3}^{4} \gamma_{\epsilon}(l) + \sum_{l=3}^{6} \gamma_{\epsilon}(l)) \sum_{i=1}^{2} \gamma_{\epsilon}(i)}, $$

$$ \phi^{\pm, (8)} = \phi^{\pm, (0)} |\mathcal{I}^{\epsilon}(t; a) \phi_0|^{\sum_{i=1}^{8} \gamma_{\epsilon}(i) + \sum_{l=7}^{8} \gamma_{\epsilon}(l) \sum_{i=1}^{6} \gamma_{\epsilon}(i) + (\sum_{l=5}^{6} \gamma_{\epsilon}(l) + \sum_{l=5}^{8} \gamma_{\epsilon}(l)) \sum_{i=1}^{4} \gamma_{\epsilon}(i)} \\ \times |\mathcal{I}^{\epsilon}(t; a) \phi_0|^{(\sum_{l=3}^{4} \gamma_{\epsilon}(l) + \sum_{l=3}^{6} \gamma_{\epsilon}(l) + \sum_{l=3}^{8} \gamma_{\epsilon}(l) + \sum_{l=3}^{4} \gamma_{\epsilon}(l) \sum_{l=7}^{8} \gamma_{\epsilon}(l)) \sum_{i=1}^{2} \gamma_{\epsilon}(i)}. $$

By induction, we get at $(0, T) \times \Gamma_\eta^\pm$ for all $k$ large enough

$$ \phi^{\pm, (2k)} = \mathcal{I}^{\epsilon}(t; a) \phi_0 |\mathcal{I}^{\epsilon}(t; a) \phi_0|^{\sum_{i=1}^{2k} \gamma_{\epsilon}(i) + R_k}, \quad (4.13) $$

which is finite for any $k$ by $\sum_{k \in \mathbb{N}^*} \gamma_\epsilon(k)$ being assumed finite, and where we have defined

$$ R_k := \sum_{p=1}^{k-1} (\Lambda_p(k) \sum_{i=1}^{2(k-p)} \gamma_\epsilon(i)), $$

where $\{\Lambda_p\}_{p \ge 1}(k)$ is a bounded positive real function that can be explicitly determined and tends to zero faster than $\gamma_\epsilon(k)$ when $k$ goes to infinity. In particular, the first three terms are given by

$$ \left\{ \begin{array}{l} \Lambda_1(k) = \gamma_{\epsilon}(2k)\gamma_{\epsilon}(2k-1), \\ \Lambda_2(k) = \sum_{l=2k-3}^{2k-2} \gamma_{\epsilon}(l) + \sum_{l=2k-3}^{2k} \gamma_{\epsilon}(l), \\ \Lambda_3(k) = \sum_{l=2k-5}^{2k-4} \gamma_{\epsilon}(l) + \sum_{l=2k-5}^{2k-2} \gamma_{\epsilon}(l) + \sum_{l=2k-5}^{2k} \gamma_{\epsilon}(l) + \sum_{l=2k-5}^{2k-4} \gamma_{\epsilon}(l) \sum_{l=2k-1}^{2k} \gamma_{\epsilon}(l). \end{array} \right. $$
Denoting again $e^{\pm,(k)} := \phi^{\pm,(k)} - \Psi_e$, we get at $(0, T) \times \Gamma_{\eta}^{\pm}$,

$$e^{\pm,(2k)} = \mathcal{I}^{\epsilon}(t; a)\phi_0 |\mathcal{I}^{\epsilon}(t; a)\phi_0|^{\sum_{i=1}^{2k} \gamma_{\epsilon}(i) + \sum_{l=2k-1}^{2k} \gamma_{\epsilon}(l)} \sum_{i=1}^{2k-2} \gamma_{\epsilon}(i) + (\sum_{l=2k-3}^{2k-2} \gamma_{\epsilon}(l) + \sum_{l=2k-3}^{2k} \gamma_{\epsilon}(l)) \sum_{i=1}^{2k-4} \gamma_{\epsilon}(i) + R_k - \Psi_e. \quad (4.14)$$

As $\sum_i \gamma_\epsilon(i)$ is convergent, for any $\chi \in (0, 1)$, there exists $k_\chi \in \mathbb{N}^*$ such that for all $k \ge k_\chi$, we get from (4.13) at $(0, T) \times \Gamma_\eta^\pm$,

$$1 \le \frac{\phi^{\pm,(2k+2)}}{\phi^{\pm,(k)}} \le |\mathcal{I}^{\epsilon}(t; a)\phi_0|^{\chi}.$$

As a consequence, for any $\zeta > 0$ as small as wanted, there exists $k_\zeta \ge 2$, such that $|e^{\pm,(k)} - e^{\pm,(k-2)}| < \zeta$, for any $k \ge k_\zeta$. We deduce that (4.7) is asymptotically convergent as CSWR. Moreover, for $\gamma_\epsilon$ satisfying Assumption iii), and for $k_0$ as defined in Case 2, we have, in particular, for all $k \le k_0/2$ at $(0, T) \times \Gamma_\eta^\pm$

$$|e^{\pm,(2k)}| \le |c_\eta(t, \cdot)| = |\mathcal{I}^{\epsilon}(t; a)\phi_0 - \Psi_e(t, \cdot)|. \quad (4.15)$$

Then, from (4.15) and from the existence of a $c(t)$ such that $\|\mathcal{I}^{\epsilon}(t; a)\phi_0 - \Psi_e(t, \cdot)\|_{L^2(\mathbb{R}^d)} \le c(t)\epsilon\|\phi_0\|_{L^2(\mathbb{R}^d)}$ [39], then

$$\|\phi^{(2k)}(t, \cdot) - \Psi_e(t, \cdot)\|_{L^2(\mathbb{R}^d)} \le \|\mathcal{I}^{\epsilon}(t; a)\phi_0 - \Psi_e(t, \cdot)\|_{L^2(\mathbb{R}^d)} \le c(t)\epsilon\|\phi_0\|_{L^2(\mathbb{R}^d)}, \quad (4.16)$$

where

$$\phi^{(k)}(t, \mathbf{y}) = \begin{cases} \phi^{\pm,(k)}(t, \mathbf{y}), & \text{for } (t, \mathbf{y}) \in (0, T) \times \Omega_{\eta}^{\pm} \setminus (\Omega_{\eta}^{+} \cap \Omega_{\eta}^{-}), \\ \frac{\phi^{+, (k)} + \phi^{-, (k)}}{2}(t, \mathbf{y}), & \text{for } (t, \mathbf{y}) \in (0, T) \times \Omega_{\eta}^{+} \cap \Omega_{\eta}^{-}. \end{cases}$$

Due to (4.16), any given convergence tolerance will be then reached in less iterations than CSWR and than in Case 2. $\square$

### 4.5. Computational complexity

As discussed in Section 2, parallel computation of FGA can be efficiently implemented in the semi-classical regime by decomposing the summation in the phase space. This feature can be used for HKSWR in the middle frequency regime to speed-up the computation of FGA in HKSWR. It requires two levels of parallelization, one in space (interior scheme) and one in frequency (transmission conditions). As shown in Section 2, the computational complexity of FGA at a given point $\mathbf{x}_j$ is given by (2.9), denoted as $\text{Op}(\phi_{\text{FGA}}^{\epsilon}(t^n, \mathbf{x}_j))$; this can be reduced by a half by parallel implementation on two equal-sized subdomains. Let us denote by $D_{\mathbf{x}}^{\pm}$ two real space domains in HKSWR, with $N_{\mathbf{x}}^{\pm}$ degrees of freedom. Each time iteration requires about $\mathcal{O}((N_{\mathbf{x}}^{\pm})^{\alpha^{\pm}})$ operations, for some $\alpha^{\pm} > 1$. Assuming that the overlapping region $\Gamma_{\mathbf{x}}$ involves $M_{\mathbf{x}}$ nodes, a parallel implementation of FGA on two subdomains requires $\sum_{\mathbf{x}_j \in \Gamma_{\mathbf{x}}} \text{Op}(\phi_{\text{FGA}}^{\epsilon}(t^n, \mathbf{x}_j))/2$. As a consequence, unlike OSWR, the computation of the transmission conditions in HKSWR is scalable and does not require the storage of the solution at the domain boundary at all time. In other words, for fixed $M_{\mathbf{x}}$, the nonlocal transmission conditions in OSWR requires a number of operations independent of the number of processors $p$ unless elaborated parallelism is developed, while the number of operations for computing the transmission conditions in HKSWR is trivially perfectly scalable. The consequence is that the computation of FGA-based transmission condition decreases linearly with the number of processors.
## 5. Geometric Optics Schwarz Waveform Relaxation (GOSWR)

In this section, we develop another semiclassical Schwarz Waveform Relaxation method based on geometric optics, namely the Geometric Optics Schwarz Waveform Relaxation (GOSWR) method. The motivation is to further improve the computational complexity of HKSWR by replacing FGA with a rougher but cheaper approximation to $\psi^\epsilon$ in (1.1) in the middle frequency regime. Therefore the tradeoff here is that, GOSWR will accelerate computation of each Schwarz iteration but decrease the convergence rate, as tested in numerical examples in Section 6.

Recall that an N-th order geometric optics approximation (GOA) is given by the following WKB ansatz

$$ \phi_{\text{GOA}}^{\epsilon,N,\pm} = \sum_{l=0}^{N} a_l(t, \mathbf{x}) e^{iS(t, \mathbf{x})/\epsilon}. \quad (5.1) $$

Then the phase function $S$ satisfies the eikonal equation and the leading order term $a_0$ is solution to the transport equation

$$ \begin{aligned} \partial_t S + \frac{1}{2} |\nabla S|^2 + V &= 0, \\ \partial_t a_0 + \nabla S \cdot \nabla a_0 + \frac{a_0}{2} \Delta S &= 0. \end{aligned} $$

Higher order terms $a_l$ ($l \le N$) can iteratively be computed by

$$ \partial_t a_n + \nabla S \cdot \nabla a_n + \frac{1}{2} a_n \Delta S = \Delta a_{n-1}. $$

Note that (5.1) is only valid before the formation of caustics. However, from a practical point of view, we do not expect any difficulty since GOA will only be applied on a sufficiently short time (typically $\Delta t$) in Schwarz iterations. Now let us introduce GOSWR as follows

$$ \left\{ \begin{array}{ll} P_{\epsilon} \phi^{\pm,(k)}(0, \cdot) = 0, & \text{on } (0, T) \times \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(0, \cdot) = \phi_0(\cdot), & \text{on } \Omega_{\eta}^{\pm}, \\ \phi^{\pm,(k)}(t, \cdot) = \mathcal{F}_{\gamma_{\epsilon}(k)}(\phi_{\text{GOA}}^{\epsilon,N,\pm}(t, \cdot), \phi^{(k-1)}(t, \cdot)), & \text{on } (0, T) \times \Gamma_{\eta}^{\pm}, \end{array} \right. \quad (5.2) $$

which is the same as (4.11), except FGA being replaced by GOA. Specifically, one can construct GOA by the method of characteristics. For a non-zero potential, a numerical integrator is needed for computing $a^\epsilon$ and $S$. In addition, the Schwarz waveform method proposed in (5.2) only requires the knowledge of the amplitude but not the phase, and thus further reduce the computational complexity of the method. We denote by $\mathbf{q}^n = (q_j^n)_j, \mathbf{p}^n = (p_j^n)_j, S^n = (S_j^n)_j, \mathbf{a}^n = (a_j^n)_j$ respectively the particle positions, momenta, phases and amplitudes at time $t_n$ which are assumed to be known. Then from time $t_n$ to $t_{n+1}$,

1. We solve for $t \in (t_n, t_{n+1})$,

$$ \frac{d\mathbf{q}}{dt} = \mathbf{p}, \quad \frac{da}{dt} = -\frac{1}{2}a\Delta S, \quad \frac{dS}{dt} = \frac{\mathbf{p}^2}{2} - V, \quad \frac{d\mathbf{p}}{dt} = -\nabla V, $$

to get $\mathbf{a}^{n+1}$, $S^{n+1}$, $\mathbf{p}^{n+1}$ and $\mathbf{q}^{n+1}$. The use of a short time step will avoid the issue of caustics.

2. We construct $|\phi_{\text{GOA}}^{\epsilon,N,\pm}(t_{n+1}, \cdot)|^{\gamma_{\epsilon}(k)}$ at $\Gamma_{\eta}^{\pm}$ from $\mathbf{a}^{n+1}$ by (5.1).

## 6. Numerical tests

In this section, we present the performance of the semiclassical Schwarz Waveform Relaxation methods introduced in previous sections.
## 6.1. High-frequency regime

We first validate the domain decomposition approach presented in Section 2 in the high-frequency regime. We consider a one-dimensional problem ($d = 1$). Tests in higher dimension will actually be essentially the same since the parallelism structure is similar to the one-dimensional case. We choose the following initial data and potential

$$ \phi_0(x) = \exp(-25x^2 + i(x^2 + 3x)/\epsilon), \quad V(x) = \exp(-0.2x^2). $$

We determine $\phi_{FGA}$ at time $T = 1$, for $\Delta t = 10^{-5}$, $\epsilon = 2^{-9}$ and the spatial domain of computation ($-40, 40$). The FGA is computed following the procedure described above with RK4 scheme for simulating the Hamiltonian flow. Details can also be found in [43]. We represent in Figure 6.1, $\phi_{FGA}$, $\phi_{FGA}^{+}$, $\phi_{FGA}^{-}$ at final time $T = 1$, corresponding respectively to the total FGA, the contribution from the first subdomain in $\mathbf{q}$ ($\Omega_{\mathbf{q}}^{+}$ and $\Omega_{\mathbf{q}}^{-}$) and which are respectively defined in (2.8) and (2.10). Naturally, we have $\phi_{FGA} = \phi_{FGA}^{+} + \phi_{FGA}^{-}$. The complexity to compute $\phi_{FGA}$ on one domain is as expected roughly twice the one for computing $\phi_{FGA}^{\pm}$ on two disjoint subdomains.

<figure><img src="image_5.png" /><figcaption>Figure 6.1: FGA modulus computation by using two subdomains: $\phi_{FGA}^{\epsilon, \pm}$ and reconstructed global solution $\phi_{FGA}^{\epsilon} = \phi_{FGA}^{\epsilon,+} + \phi_{FGA}^{\epsilon,-}$.</figcaption></figure>

The computational times are summarized in Table 6.1, and confirm that parallel computing for FGA can be performed with an ideal speed-up.

<table><caption>Table 6.1: CPU time comparison</caption><thead><tr><td>FGA domain</td><td>CPU time (in sec.)</td></tr></thead><tbody><tr><td>Contribution Domain 1 ($\phi_{FGA}^{\epsilon,+}$)</td><td>34.1</td></tr><tr><td>Contribution Domain 2 ($\phi_{FGA}^{\epsilon,-}$)</td><td>35.2</td></tr><tr><td>One unique domain ($\phi_{FGA}^{\epsilon}$)</td><td>70.4</td></tr></tbody></table>

## 6.2. Middle-frequency regime

This subsection is devoted to numerical tests for the middle-frequency regime as described in Section 4. In the following, we compare the CSWR (4.8), OSWR (3.8), HKSWR (4.11) and GOSWR methods (5.2). For solving the particle and phase evolutions, RK4 is used in HKSWR and Heun's scheme in GOSWR.
### 6.2.1. Test I

We consider the one-dimensional Schrödinger equation ($d = 1$)

$$i\epsilon\frac{\partial\psi^{\epsilon}}{\partial t} = -\frac{\epsilon^2}{2}\Delta\psi^{\epsilon}(t, x) + V(x)\psi^{\epsilon}(t, x), x \in \Omega = (-2, 2), t \in (0, T),$$

with $\epsilon = 2^{-3}$. We choose the initial data and potential to

$$\psi_0(x) = \exp(-25x^2 + i(x^2 + 3x)/\epsilon), \quad V(x) = 10\exp(-0.2x^2).$$

The Schrödinger equation is discretized with a second-order finite difference scheme on $\Omega_{\eta}^{\pm}$. The parameters for the simulation are as follow. Respective space and time steps are $\Delta x = 2^{-5}$, $\Delta t = 5 \times 10^{-3}$. The final physical time is $T = 0.1$. The overlapping region is reduced to $\eta = \Delta x$, and more specifically $\Omega_{\eta}^{+} = (-2, \Delta x)$ and $\Omega_{\eta}^{-} = (0, 2)$. The numerical simulations are performed by choosing $\tilde{\epsilon} = 2^{-6}$ and $\gamma_{\epsilon}(k) = \exp(-k^2/20)$, where $k$ refers to as the Schwarz iteration.

In Figure 6.2, we report, with respect to the Schwarz iterations, i) the $\ell^2$-norm in time of $\phi^{+}(t, \Delta x) - \phi^{-}(t, \Delta x)$ and ii) the $\ell^2$-error at time $T = 0.1$ in space, between the reconstructed solution and a solution of reference. We compare the convergence rate of HKSWR (4.11) with the simple CSWR algorithm (4.8) as well as OSWR as described in [12], based on second-order pseudodifferential absorbing boundary conditions. More specifically: "Order 2 OSWR" refers to as transmission conditions based on a DtN operator including effects of the potential, which writes at $(0, T) \times \Gamma_{\eta}^{\pm}$,

$$(\pm \sqrt{\epsilon}\partial_x + B_{\epsilon}^{(1)})\phi^{\pm,(k)} = (\pm \sqrt{\epsilon}\partial_x + B_{\epsilon}^{(1)})\phi^{\mp,(k-1)},$$

where $B_{\epsilon}^{(1)} := e^{i\pi/4}\partial_t^{1/2} - tV'(x) - iV'(x)I_t/4$. The operator $\partial_t^{1/2}$ is the half-order derivative operator while $I_t$ denotes the integral operator. Details can be found in [12].

<figure><img src="image_6.png" /><figcaption>Figure 6.2: i) $\ell^2$-error in time at the boundary (first interface) $\|\phi^{+}(0, \cdot) - \phi^{-}(0, \cdot)\|_{L^2(0,T)}$ as a function of Schwarz iterations, ii) $\ell^2$-error in space, as a function of Schwarz iteration at time $T = 0.1$ and $\epsilon = 2^{-3}$.</figcaption></figure>

If $\epsilon$ is not small enough as in this example, it is not surprising to observe that transmission conditions based on TBCs or high-order ABCs allow for a faster convergence than HKSWR and GOSWR. Indeed, the FGA does not provide a very accurate approximation of the Schrödinger equation for $\epsilon = 2^{-3}$. As a consequence, although it is improved compared to CSWR, the HKSWR and GOSWR algorithms still have a slower convergence rate than OSWR. The expectation is however that for a smaller values of $\epsilon$, (4.11) and (5.2) approaches will be more relevant from an efficiency point of view.
### 6.2.2. Test II

This test is similar to the one presented in the numerical test I, except that $\epsilon$ is now taken smaller: $\epsilon = 2^{-5}$, (resp. $\epsilon = 2^{-6}$). The Schrödinger equation is discretized with finite differences on $\Omega_{\eta}^{\pm}$ for the following simulation data. Respective space and time steps are $\Delta x = 2^{-7}$ (resp. $\Delta x = 2^{-9}$), $\Delta t = 5 \times 10^{-4}$ (resp. $\Delta t = 10^{-4}$). The final time of computation is $T = 0.1$. The overlapping region is reduced to $\eta = \Delta x$, and more specifically $\Omega_{\eta}^{-} = (-2, \Delta x)$ and $\Omega_{\eta}^{+} = (0, 2)$. Numerical tests are performed by choosing $\tilde{\epsilon} = 2^{-10}$ and $\gamma_{\epsilon}(k) = \exp(-(k-1)^2/100)$, where $k$ denotes the Schwarz iteration. Figure 6.3 reports, with respect to the Schwarz iteration $k$, (resp. Figure 6.5) i) the $l^2$-norm in time of $\phi^{+}(t, \Delta x) - \phi^{-}(t, \Delta x)$, and ii) the $l^2$-error in space of the reconstructed solution with the solution of reference computed at final time $T = 0.1$.

<figure><img src="image_7.png" /><figcaption>Figure 6.3: i) $l^2$-error in time at the boundary (first interface) $\|\phi^{+}(0, \cdot) - \phi^{-}(0, \cdot)\|_{L^2(0,T)}$ as a function of Schwarz iterations, ii) $l^2$-error in space, as a function of Schwarz iteration at time $T = 0.1$ and $\epsilon = 2^{-5}$.</figcaption></figure>

We also represent the reference and reconstructed solutions on Figure 6.4, for $\epsilon = 2^{-7}$.

<figure><img src="image_8.png" /><figcaption>Figure 6.4: Comparison of the i) modulus, ii) real part, iii) imaginary part of the reconstructed solution $(\phi^{-}, \phi^{+})$ with the amplitude of the reference solution (computed in all $\Omega$) at time $T = 0.1$. The domain interface is located at $x = 0$ and $\epsilon = 2^{-7}$.</figcaption></figure>

We observe that the smaller $\epsilon$, the slower the convergence of the FGA-like DDM, Figures 6.2, 6.3, 6.5. These results also suggest that GOA is a good compromise between computational complexity and rate of convergence. As expected, the convergence rate is indeed not as fast as FGA, but it has a lower computational complexity and is faster than CSWR.
<figure><img src="image_9.png" /><figcaption>Figure 6.5: i) $ℓ^2$-error in time at the boundary (first interface) $||\phi^+(0, \cdot) - \phi^-(0, \cdot)||_{L^2(0,T)}$ as a function of Schwarz iterations, ii) $ℓ^2$-error in space, as a function of Schwarz iteration at time $T = 0.1$ and $\varepsilon = 2^{-6}$.</figcaption></figure>

## 7. Conclusion

This paper was devoted to the derivation of domain decomposition methods for the time-dependent Schrödinger equation in the middle/high-frequency regime by using Schwarz waveform relaxation methods with FGA-based Dirichlet transmission conditions. The idea is to simultaneously benefit from the simplicity of the CSWR algorithm and the accuracy of FGA. Convergence of the so-called Herman-Kluk Schwarz Waveform Relaxation algorithm (HKSWR) is indeed a direct consequence of the CSWR convergence, but with acceleration provided by FGA. Some numerical evidences were proposed showing the efficiency of the method. Table 7.1 roughly summarizes and compares the behavior of each presented method depending on the regime. Among the criteria from Table 7.1, *MultiD easyness* refers to the easyness to implement, in particular in parallel, the corresponding DDM in more than one dimension. Although this criterion is quite subjective, we can surely ensure that the complexity of implementation of FGA-based DDM (FGA, FGA-based, HKSWR) or CSWR is independent of the spatial dimension, which is not true for DtN-based transmission conditions (OSWR). *Comput. Comp.* refers to as the computational cost to implement the transmission condition, by comparison with simple Dirichlet-based transmission conditions (CSWR).

<table><caption>Table 7.1: DDM comparison</caption><thead><tr><th></th><th>FGA</th><th>FGA-based</th><th>CSWR</th><th>OSWR</th><th>HKSWR</th><th>GOSWR</th></tr><tr><th></th><th>Algo. (2.10)</th><th>Algo. (4.4)</th><th>Algo. (4.8)</th><th>Algo. (3.8)</th><th>Algo. (4.11)</th><th>Algo. (5.2)</th></tr></thead><tbody><tr><th>Semi-classical Regime</th><td>Y</td><td>N</td><td>N</td><td>N</td><td>N</td><td>N</td></tr><tr><th>Quantum Regime</th><td>N</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><th>Convergence</th><td>Y</td><td>N</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><th>Rate of Convergence</th><td>Fast</td><td>-</td><td>Slow</td><td>Fast</td><td>Fast</td><td>Fast</td></tr><tr><th>Comput. Comp.</th><td>Y</td><td>Y</td><td>N</td><td>Y</td><td>Y</td><td>N</td></tr><tr><th>MultiD Easyness</th><td>Y</td><td>Y</td><td>Y</td><td>N</td><td>Y</td><td>Y</td></tr></tbody></table>

The HKSWR and GOSWR algorithms are planed to be implemented in higher dimension and tested on more realistic problems, in particular related to laser-molecule interactions which involve high frequency generation by electronic recombination [12, 14, 31]. We conclude by mentioning that the HKSWR method presented here can in principle be extended to hyperbolic systems of conservation laws and wave equations using the corresponding FGA [34, 36].
Acknowledgements.

E. Lorin authors thanks NSERC for the financial support via the Discovery Grant program. X. Yang was partially supported by the NSF grants DMS-1418936, and DMS-1107291: NSF Research Network in Mathematical Sciences “Kinetic description of emerging challenges in multiscale problems of natural science”, and Hellman Family Foundation Faculty Fellowship, UC Santa Barbara. X. Antoine thanks the support of the french ANR grants ”Bond” (ANR-13-BS01-0009-01).

References

[1] X. Antoine, A. Arnold, C. Besse, M. Ehrhardt, and A. Schaedle. A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations. *Communications in Computational Physics*, 4(4):729–796, 2008.

[2] X. Antoine and H. Barucq. Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism. *SIAM J. Appl. Math.*, 61(6):1877–1905 (electronic), 2001.

[3] X. Antoine and H. Barucq. On the construction of approximate boundary conditions for solving the interior problem of the acoustic scattering transmission problem. In *Domain decomposition methods in science and engineering*, volume 40 of *Lect. Notes Comput. Sci. Eng.*, pages 133–140. Springer, Berlin, 2005.

[4] X. Antoine and C. Besse. Construction, structure and asymptotic approximations of a microdifferential transparent boundary condition for the linear Schrödinger equation. *J. Math. Pures Appl. (9)*, 80(7):701–738, 2001.

[5] X. Antoine and C. Besse. Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation. *J. Comput. Phys.*, 188(1):157–175, 2003.

[6] X. Antoine, C. Besse, and S. Descombes. Artificial boundary conditions for one-dimensional cubic nonlinear Schrödinger equations. *SIAM J. Numer. Anal.*, 43(6):2272–2293 (electronic), 2006.

[7] X. Antoine, C. Besse, and P. Klein. Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential. *J. Comput. Phys.*, 228(2):312–335, 2009.

[8] X. Antoine, C. Besse, and P. Klein. Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential. Part I: Construction and a priori estimates. *Math. Models Methods Appl. Sci.*, 22(10):1250026, 38, 2012.

[9] X. Antoine, C. Besse, and P. Klein. Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential. Part II: Discretization and Numerical Results. *Numerische Mathematik*, 125(2):191–223, 2013.

[10] X. Antoine, C. Besse, and J. Szeftel. Towards accurate artificial boundary conditions for nonlinear PDEs through examples. *Cubo*, 11(4):29–48, 2009.

[11] X. Antoine and E. Lorin. Asymptotic estimates of the convergence of schwarz waveform relaxation domain decomposition methods for stationary classical and quantum waves. *To Be Submitted*, 2015.

[12] X. Antoine, E. Lorin, and A.D. Bandrauk. Domain decomposition method and high-order absorbing boundary conditions for the numerical simulation of the time dependent Schrödinger equation with ionization and recombination by intense electric field. *Journal of Scientific Computing*, 2014.

[13] X. Antoine, E. Lorin, J. Sater, F. Fillion-Gourdeau, and A. D. Bandrauk. Absorbing boundary conditions for relativistic quantum mechanics equations. *J. Comput. Phys.*, 277:268–304, 2014.
[14] A. D. Bandrauk, F. Fillion-Gourdeau, and E. Lorin. Atoms and molecules in intense laser fields: gauge invariance of theory and models. *Journal of Physics B-Atomic Molecular and Optical Physics*, 46(15), AUG 14 2013.

[15] W. Bao, S. Jin, and P.A. Markowich. On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime. *J. Comput. Phys.*, 175:487-524, 2002.

[16] J.-P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. *Journal of Computational Physics*, 114(2):185-200, 1994. cited By 5129.

[17] V. Dolean, M. J. Gander, and L. Gerardo-Giorda. Optimized Schwarz methods for Maxwell's equations. *SIAM J. Sci. Comput.*, 31(3):2193-2213, 2009.

[18] M. El Bouajaji, V. Dolean, M. J. Gander, and S. Lanteri. Optimized Schwarz methods for the time-harmonic Maxwell equations with damping. *SIAM J. Sci. Comput.*, 34(4):A2048-A2071, 2012.

[19] B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. *Math. Comp.*, 31(139):629-651, 1977.

[20] B. Engquist and O. Runborg. Computational high frequency wave propagation. *Acta Numer.*, 12:181-266, 2003.

[21] M. Gander and L. Halpern. Optimized schwarz waveform relaxation methods for advection reaction diffusion problems. *SIAM J. Num. Anal.*, 45(2), 2007.

[22] L. Halpern and J. Rauch. Error analysis for absorbing boundary conditions. *Numer. Math.*, 51(4):459-467, 1987.

[23] L. Halpern and J. Rauch. Absorbing boundary conditions for diffusion equations. *Numer. Math.*, 71(2):185-224, 1995.

[24] L. Halpern and J. Szeftel. Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation. *Math. Models Methods Appl. Sci.*, 20(12):2167-2199, 2010.

[25] M.F. Herman and E. Kluk. A semiclassical justification for the use of non-spreading wavepackets in dynamics calculations. *Chemical Physics*, 91(1):27-34, 1984. cited By (since 1996)484.

[26] L. Hörmander. *Linear partial differential operators*. Springer Verlag, Berlin, 1976.

[27] L. Hörmander. *The analysis of linear partial differential operators. III*. Classics in Mathematics. Springer, Berlin, 2007. Pseudo-differential operators.

[28] S. Jin, P.A. Markowich, and C. Sparber. Mathematical and computational methods for semiclassical Schrödinger equations. *Acta Numer.*, 20:211-289, 2011.

[29] S. Jin, H. Wu, and X. Yang. Gaussian beam methods for the Schrödinger equation in the semi-classical regime: Lagrangian and Eulerian formulations. *Commun. Math. Sci.*, 6:995-1020, 2008.

[30] R. Kechroud, X. Antoine, and A. Soulaïmani. Numerical accuracy of a Padé-type non-reflecting boundary condition for the finite element solution of acoustic scattering problems at high-frequency. *Internat. J. Numer. Methods Engrg.*, 64(10):1275-1302, 2005.

[31] M. Lewenstein, Ph. Balcou, M.Yu. Ivanov, A. L'Huillier, and P.B. Corkum. Theory of high-harmonic generation by low-frequency laser fields. *Phys. Rev. A*, 49(3):2117-2132, 1994.

[32] J. Li, G. Lin, and X. Yang. A frozen gaussian approximation-based multi-level particle swarm optimization for seismic inversion, 2015. in revision.
[33] E. Lorin, S. Chelkowski, and A. D. Bandrauk. Mathematical modeling of boundary conditions for laser-molecule time-dependent Schrödinger equations and some aspects of their numerical computation—one-dimensional case. *Numer. Methods Partial Differential Equations*, 25(1):110–136, 2009.

[34] J. Lu and X. Yang. Frozen Gaussian approximation for high frequency wave propagation. *Commun. Math. Sci.*, 9(3):663–683, 2011.

[35] J. Lu and X. Yang. Convergence of frozen Gaussian approximation for high-frequency wave propagation. *Comm. Pure Appl. Math.*, 65(6):759–789, 2012.

[36] J. Lu and X. Yang. Frozen Gaussian approximation for general linear strictly hyperbolic systems: formulation and Eulerian methods. *Multiscale Model. Simul.*, 10(2):451–472, 2012.

[37] P. Markowich, P. Pietra, and C. Pohl. Numerical approximation of quadratic observable of Schrödinger equation-type equations in the semiclassical limit. *Numer. Math.*, 81:595–630, 1999.

[38] P. Markowich, P. Pietra, C. Pohl, and H. Stimming. A Wigner-measure analysis of the Dufort-Frankel scheme for the Schrödinger equation. *SIAM J. Numer. Anal.*, 40:1281–1310, 2000.

[39] T. Swart and V. Rousse. A mathematical justification for the Herman-Kluk propagator. *Comm. Math. Phys.*, 286(2):725–750, 2009.

[40] J. Széftel. Absorbing boundary conditions for one-dimensional nonlinear Schrödinger equations. *Numer. Math.*, 104(1):103–127, 2006.

[41] M. E. Taylor. *Partial differential equations III. Nonlinear equations*, volume 117 of *Applied Mathematical Sciences*. Springer, New York, second edition, 2011.

[42] V. Vaibhav. Artificial boundary conditions for certain evolution PDEs with cubic nonlinearity for non-compactly supported initial data. *J. Comput. Phys.*, 230(8):3205–3229, 2011.

[43] X. Yang and J. Zhang. Computation of the Schrödinger equation in the semiclassical regime on an unbounded domain. *SIAM J. Numer. Anal.*, 52(2):808–831, 2014.