X Antoine
email: xavier.antoine@univ-lorraine.fr

E Lorin
email: elorin@math.carleton.ca

On the rate of convergence of Schwarz waveform relaxation methods for the time-dependent Schrödinger equation

Keywords: Schrödinger equation, domain decomposition method, pseudo-differential calculus

This paper is dedicated to the analysis of the rate of convergence of the classical and quasioptimal Schwarz waveform relaxation (SWR) method for solving the linear Schrödinger equation with space-dependent potential. The strategy is based on i) the rewriting of the SWR algorithm as a fixed point algorithm in frequency space, and ii) the explicit construction of contraction factors thanks to pseudo-differential calculus. Some numerical experiments illustrating the analysis are also provided.

Introduction

This paper is devoted to the analysis of the rate of convergence of a class of domain decomposition methods (DDM), the Schwarz Waveform Relaxation algorithms for solving the real-time linear Schrödinger equation (LSE) [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]. The analyzed DDMs provide algorithms for solving all kinds of wave propagation problems, and mainly involve two concepts: relaxation algorithms and well-designed transmission conditions [START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF].

Let us consider the following initial boundary-value problem: determine the complexvalued wavefunction u solution to the LSE

i∂ t u + u + V (x)u = 0, x ∈ R, t ≥ 0, |u(x, t)| → x→±∞ 0, t ≥ 0, u(x, 0) = u 0 (x), x ∈ R, (1)
where u 0 is the Cauchy data, V is a real-valued space-dependent smooth potential, which is positive (respectively negative) for attractive (respectively repulsive) interactions. A Schwarz Waveform Relaxation (SWR) domain decomposition algorithm is proposed to solve the initial value problem [START_REF] Al-Khaleel | Optimized waveform relaxation methods for longitudinal partitioning of transmission lines[END_REF]. In the one-dimensional case, the SWR algorithm decomposes R into two (or generally more) regions, with or without overlap, and two uncoupled IBVPs are then considered. The SWR iterations allow for a reconstruction of the solution between the adjacent subdomains by exchanging information at the subdomain interfaces thanks to the transmission conditions. This makes these IBVs on each subdomain suitable for parallel computing [START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation[END_REF].

In this paper, we analyze the rate of convergence of the Classical and quasi-Optimal Schwarz Waveform Relaxation (CSWR and q-OSWR) DDMs [START_REF] Al-Khaleel | Optimized waveform relaxation methods for longitudinal partitioning of transmission lines[END_REF][START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation[END_REF][START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Gander | Optimal convergence for overlapping and nonoverlapping Schwarz waveform relaxation[END_REF][START_REF] Gander | Optimal Schwarz waveform relaxation methods for the one-dimensional wave equation[END_REF][START_REF] Gander | Dirichlet-neumann and neumann-neumann waveform relaxation for thewave equation[END_REF][START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF][START_REF]A time-dependent Dirichlet-Neumann method for the heat equation[END_REF], applied to the time-dependent Schrödinger equation in one-dimension. Although these methods have received much attention over the past decades, the first application to the Schrödinger equation can be found in [START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF], where the authors consider the real-time linear one-dimensional Schrödinger equation. In another recent paper [START_REF] Antoine | Domain decomposition method and highorder absorbing boundary conditions for the numerical simulation of the time dependent Schrödinger equation with ionization and recombination by intense electric field[END_REF], some algorithms are analyzed for the one-dimensional time-dependent Linear Schrödinger Equation (LSE) where are included ionization and recombination processes by an intense electric field, and in [START_REF] Lorin | Schwarz waveform relaxation domain decomposition methodology for the N-body time-independent and time-dependent Schrödinger equation[END_REF] a SWR methodology for solving the N -body Schrödinger equation is considered. In [START_REF] Besse | Schwarz waveform relaxation method for one dimensional Schrödinger equation with general potential[END_REF], the authors study the numerical performance with a GPU implementation of Schwarz waveform relaxation methods for the one-dimensional dynamical solution of the LSE with a general potential. It is shown that these algorithms are fast and robust for complex linear problems. In [START_REF] Lorin | Frozen gaussian approximation based domain decomposition methods for the linear and nonlinear Schrödinger equation beyond the semi-classical regime[END_REF], domain decomposition methods have been developed and combined with geometric optics and frozen gaussian approximation for computing the solution to the linear Schrödinger equations under and beyond the semi-classical regime. Finally, [START_REF] Antoine | Lagrange -Schwarz waveform relaxation domain decomposition methods for linear and nonlinear quantum wave problems[END_REF] is dedicated to the development of high-order transmission conditions for SWR methods applied to the Schrödinger equation, using only local operators. The goal of the present paper is to contribute to the mathematical understanding of Schwarz waveform relaxation DDMs for solving the time-dependent Schrödinger equation. Although SWR-DDM methods are now extensively used in all kinds of high dimensional problems, the rigorous analysis of the rate of convergence remains incomplete, and is only understood in some simple configurations. This paper is precisely dedicated to this question for the Schrödinger equation in the one-dimensional case with non-constant potentials. The strategy which is proposed can in principle, be applied in higher dimension and in the stationary case, see [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF][START_REF] Antoine | Asymptotic estimates of the convergence of classical Schwarz waveform relaxation domain decomposition methods for two-dimensional stationary quantum waves[END_REF]. It extensively uses pseudo-differential calculus [START_REF] Alinhac | Pseudo-differential operators and the Nash-Moser theorem[END_REF][START_REF] Hörmander | Linear Partial Differential Operators[END_REF][START_REF] Nirenberg | Lectures on Linear Partial Differential Equations[END_REF], and was originally developed for deriving and analyzing high-order absorbing boundary conditions for classical and quantum wave equations, as well as diffusion equations [START_REF] Antoine | Absorbing boundary conditions for the twodimensional Schrödinger equation with an exterior potential. Part I: construction and a priori estimates[END_REF][START_REF] Antoine | Absorbing boundary conditions for the twodimensional Schrödinger equation with an exterior potential. Part II: discretization and numerical results[END_REF][START_REF] Antoine | A friendly review to absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum wave equations[END_REF][START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF][START_REF] Halpern | Error analysis for absorbing boundary conditions[END_REF]. Accurate transmission operators (which will provide fast convergence) can indeed be derived by using asymptotic expansions of pseudo-differential operators (in the sense of classical symbols [START_REF] Hörmander | Linear Partial Differential Operators[END_REF]). However, pseudo-differential calculus is also essential for analyzing the rate of convergence of more elementary DDMs, such as the CSWR and Robin-SWR methods which are respectively based on Dirichlet or Robin transmission operators. In Subsections 1.1 and 1.2, we introduce preliminary notations, definitions and we recall some existing results about the convergence of SWR methods. The rate of convergence of both the CSWR and the q-OSWR methods with constant and space-dependent potentials V , is analyzed in Section 2. The analysis uses pseudo-differential operator theory associated to their asymptotic symbolic calculus. In Section 3, we numerically validate the convergence rates theoretically established. We conclude in Section 4.

Brief description of the Schwarz Waveform Relaxation algorithm

In this section, we briefly describe the Classical (resp. quasi-Optimal) Schwarz Waveform Relaxation Algorithm (CSWR) (resp. (q-OSWR)). We consider a d-dimensional partial differential equation P φ = f in the spatial domain Ω ⊂ R d , and time domain (0, T ). The initial data is denoted by φ 0 . We first split Ω into two subdomains Ω ± ε with smooth boundary, with or without overlap (

Ω + ε ∩ Ω - ε = ∅ or Ω + ε ∩ Ω - ε = ∅), with > 0.
The CSWR algorithm consists in iteratively solving IBVPs in Ω ± ε × (0, T ), using Dirichlet transmission conditions at the subdomain interfaces Γ ± ε := ∂Ω ± ε , where the imposed conditions are established using the preceding Schwarz iteration data in the adjacent subdomain. For k ≥ 1, we therefore set

   P φ ±,(k) = f, in Ω ± ε × (0, T ), φ ±,(k) (•, 0) = φ ± 0 , in Ω ± ε , φ ±,(k) = φ ±,(k-1) , on Γ ± ε × (0, T ), (2)
with a given initial guess φ ±,(0) . The quasi-Optimal Schwarz Waveform Relaxation algorithm is described in the same manner as the CSWR algorithm, except that transparent or high-order absorbing transmission conditions are imposed thanks to a (pseudo-)differential operator B ± , see [START_REF] Antoine | Absorbing boundary conditions for the onedimensional Schrödinger equation with an exterior repulsive potential[END_REF]. Consequently, for k ≥ 1, we define

   P φ ±,(k) = f, in Ω ± ε × (0, T ), φ ±,(k) (•, 0) = φ ± 0 , in Ω ± ε , B (k) ± φ ±,(k) = B (k-1) ± φ ∓,(k-1) , on Γ ± ε × (0, T ), (3)
with a given initial guess φ ±,(0) .

Well-posedness of the SWR algorithm

We now define the real time-dependent one-dimensional (d = 1) Schrödinger operator P by

P (x, t, ∂ x , ∂ t ) = i∂ t + ∂ 2 x + V (x). ( 4
)
We decompose the spatial domain Ω = R into two overlapping subdomains

Ω ± ε with Ω + ε = (-∞, ε/2) and Ω - ε = (-ε/2, ∞), with ε > 0.
The Schwarz waveform relaxation with overlap corresponds to iteratively solving two IBVPs in Ω ± ε × (0, T ), by using some transmission conditions at the interfaces x = ±ε/2. The CSWR algorithm for the iteration index k = 1, 2, . . . is thus given by

   P φ ±,(k) = 0, in Ω ± ε × (0, T ), φ ±,(k) (•, 0) = φ ± 0 , in Ω ± ε , φ ±,(k) ± ε/2, • = φ ∓,(k-1) ± ε/2, • , on (0, T ), (5)
where φ ± 0 denotes the restriction of φ 0 to Ω ± ε . Introducing the error function for arbitrary iteration k and subsequently omitted, i.e. e C,± P

:= φ |Ω ± ε -φ ± , the CSWR in L 2 (Ω ± ε ) reads P e C,± P = 0 in Ω ± ε × (0, T ), e C,± P ± ε/2, t = h ± ε (t) at {±ε/2} × (0, T ), ( 6
)
where P is given by ( 4) and h ± ε are given time-dependent functions. The index P in e C,± P specifies the operator to which the error is associated to, and the upper index C stands for the CSWR algorithm (we will later use the upper index O for the q-OSWR algorithm). The following theorem states the convergence of the CSWR algorithm for a bounded potential V .

Theorem 1.1. Let ε > 0, V be in L ∞ (Ω - ε ∪ Ω + ε ). Then the algorithm (5) defines a sequence of iterates (φ +,(k) , φ -,(k) ) in H 2,1 (Ω + ε/2 × (0, T )) × H 2,1 (Ω - ε × (0, T )) with φ +,(k) (-ε/2, •), ∂ x φ +,(k) (-ε/2, •), φ -,(k) (ε/2, •) and ∂ x φ -,(k) (ε/2, •) in H 1 (0, T ).
We recall that H r,s (Ω × (0, T )) := L 2 (0, T ; H r (Ω)) ∩ H s (0, T ; L 2 (Ω)) is an anisotropic Sobolev space. This result is proven in [START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF], and is a consequence of the Trace Theorem [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]. In this paper, we are more specifically interested in estimating the convergence rate of the SWR methods. The convergence rate appears as a contraction factor, when the SWR method is rewritten as a fixed point problem. We do not recall here the details of this technical question, but we refer to [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF][START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Gander | Optimal Schwarz waveform relaxation methods for the one-dimensional wave equation[END_REF] and Subsection 2.3. Basically, in the following, we explicitly compute the contraction factor from which we can deduce the rate of convergence of the considered SWR methods.

Convergence rate of Schwarz Waveform Relaxation algorithm

Considering again the Schrödinger operator P defined in (4), we estimate in this section the theoretical convergence rates of the CSWR and q-OSWR methods. With this aim, we first provide a local factorization of the operator P in term of an incoming and outgoing wave operators.

Proposition 2.1. The Nirenberg-like factorization

P (x, t, ∂ x , ∂ t ) = (∂ x + iΛ -)(∂ x + iΛ + ) + R (7)
holds, where R ∈ OPS -∞ = m OPS m is a smooth pseudo-differential operator. The operators Λ ± are pseudo-differential operators of order 1/2 in time, and order 0 in space. In addition, their total symbol λ ± can be expanded in S 1/2 as

λ ± ∼ ∞ j=0 λ 1/2-j/2 , (8)
where λ 1/2-j/2 are elementary symbols corresponding to operators of order 1/2j/2, j ∈ N.

We can approximate λ ± by computing a finite number of inhomogeneous symbols [START_REF] Antoine | Absorbing boundary conditions for the onedimensional Schrödinger equation with an exterior repulsive potential[END_REF]. Denoting by τ the co-variable associated to t, we have: Proposition 2.2. Let us fix the principal symbol to

λ ± 1/2 = ∓ -τ + V (x) . (9)
Then, the following elementary symbols are given by

λ ± 0 = 0, λ ± -1/2 = 0 and λ ± -1 = ∓ i 4 V (x) V (x) -τ . ( 10
)
The above propositions are proved [START_REF] Antoine | Absorbing boundary conditions for the onedimensional Schrödinger equation with an exterior repulsive potential[END_REF] by using the following recursive formula, for j

∈ N * λ -j/2 = 1 2λ + 1/2 -i∂ x λ + 1/2-j/2 - j k=1 λ -j/2+k/2 λ 1/2-k/2 . ( 11
)
The following proposition which is proven in [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF], allows for a fine derivation of the transmission conditions in the q-OSWR method.

S := a ∼ +∞ j=0 a 1/2-j/2 ∈ S 1/2 S such that: a 1/2-j/2 (x, τ ) := 1 (λ + 1/2 ) j-1 L j =0 F V,1/2-j/2 (λ + 1/2 ) , with F V,1/2-j/2 ∈ C ∞ (R; R), L j ∈ N , (12)
and the associated class of pseudo-differential operators OPS 1/2 S . In [START_REF] Antoine | Domain decomposition method and highorder absorbing boundary conditions for the numerical simulation of the time dependent Schrödinger equation with ionization and recombination by intense electric field[END_REF], F V,1/2-j/2 are smooth functions depending on x and V . Then, Λ ± are in OPS

1/2
S and, for each j ∈ N, there exist some regular functions F

V,1/2-j/2 L j =0 such that λ + 1/2-j/2 = -λ - 1/2-j/2 = 1 (λ + 1/2 ) j-1 L j =0 F V,1/2-j/2 (λ + 1/2 ) . (13)

Asymptotic estimates of the contraction factor of the CSWR algorithm

In this subsection, we analyze the convergence rate of the CSWR algorithm described in [START_REF] Antoine | Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF]. As proposed in [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF], we are required to determine the contraction factor

C C P,ε of G C2 P (setting G C2 P := G C P • G C P )
, where the mapping G C P is defined from ( 6) by

G C P : h + ε , h - ε → e C,- P ε/2, • , e C,+ P -ε/2, • . (14)
The time-dependent functions h ± ε are assumed to be given. In order to simplify the notation, we denote by h ± ε the extension of h ± ε to all R which vanishes on Ω ∓ ε . We solve (6) directly to prove that G C2 P is a contraction in the (x, τ )-coordinates for V constant. For a nonconstant V , we estimate the rate of convergence through approximations. Let us then start by assuming that V is a constant (this includes V = 0). According to [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF], for a fixed time T , G C P is defined on

H 3/4 0 (0, T ) = {φ ∈ H 3/4 (0, T ) : φ(0) = 0}.
Let us characterize the part of the error e C,+ P (resp. e C,- P ) defined in [START_REF] Antoine | Absorbing boundary conditions for the onedimensional Schrödinger equation with an exterior repulsive potential[END_REF], corresponding to a traveling wave in the overlapping region Ω + ε (resp. Ω - ε ) and transmitted to the left (resp. right

) subdomain R/Ω + ε (resp. R/Ω - ε )
. Therefore, we introduce the equation

(∂ x + iΛ ∓ )e C,± Λ = 0, in Ω ± ε , e C,± Λ ± ε/2, t = h ± ε (t) at {±ε/2} × R. (15)
For V = 0, the solution to ( 15) can be explicitly computed using the Fourier transform F t along the t-direction, that is at the symbol level with respect to τ (associated to t). This exact solution to ( 15) is given in the (x, τ )-space by

êC,± Λ (x, τ ) = h ± ε (τ ) exp -i x ±ε/2 λ ∓ (τ )dy . ( 16
)
The application of Proposition (2.2) for V constant leads to:

λ ± (τ ) = λ ± 1/2 (τ ) = ∓ √ -τ + V . Then, if we define G C Λ : h + ε , h - ε → e C,- Λ ε/2, • , e C,+ Λ -ε/2, • , (17)
we have

F t G C2 Λ h + ε , h - ε = exp i ε/2 -ε/2 λ -(τ ) -λ + (τ ) dy ĥ+ ε , exp i ε/2 -ε/2 λ -(τ ) -λ + (τ ) dy ĥ- ε = exp -2iελ + (τ ) ĥ+ ε , ĥ- ε . (18)
Following [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Antoine | Absorbing boundary conditions for the onedimensional Schrödinger equation with an exterior repulsive potential[END_REF], we deduce the contraction factors C C Λ,ε of G C2 Λ in the elliptic and hyperbolic zones. The convergence of the CSWR method is ensured by the frequencies in the elliptic zones, as the frequencies in the hyperbolic zone do not affect the convergence process. In the following, we mainly restrict the analysis of the contraction factor to the elliptic zone. We denote by C C P,ε the contraction factor of G C2 P in the elliptic zone. We get:

C C P,ε = C C Λ,ε = sup τ ∈Eτ L C Λ,ε (τ ),
where

L C Λ,ε (τ ) = exp -2iελ + (τ ) = exp -2ε √ τ -V .
In the above expression the elliptic zone is denoted by

E τ = {τ ∈ R : τ > V }, and H τ is the hyperbolic zone {τ ∈ R : τ < V }.
It is also noteworthy that it is well-known that the CSWR does not converge without overlap, i.e. when ε = 0. We now consider the space-dependent potential V (x), where the hyperbolic (resp. elliptic) zone is now defined by {τ ∈ R : τ < V (x)} (resp. {τ ∈ R : τ > V (x)}). We need to characterize the error that travels from one domain to the other. For a general potential V (x), the analysis cannot be exact due to scattering effects. However, we still consider the system

(∂ x + iΛ ∓ )e C,± Λ = 0, in Ω ± ε , e C,± Λ ± ε/2, t = h ± ε (t) at {±ε/2} × R, (19)
where e C,+ Λ (resp. e C,- Λ ) is to be understood as the part of e C,+ P (resp. e C,- P ) which travels to the right (resp. left). As a consequence, the computation of e C,± Λ provides an approximation of e C,± P the solution to P e C,± P = 0. We can obtain an approximation of the contraction factor

C C P,ε , which is the contraction factor of G C2 P by C C Λ,ε for G C2 Λ C C P,ε ≈ C C Λ,ε .
This approximation will be discussed at the end of this subsection. For solving [START_REF] Gander | Optimized Schwarz methods[END_REF], let us consider the symbolic equation

(∂ x + iλ ∓ (x, τ )) e C,± Λ (x, τ ) = 0, in Ω ± ε , e C,± Λ ± ε/2, τ = h ± ε (τ ) at {±ε/2} × R. (20)
A direct computation gives

êC,± Λ (x, τ ) = h ± ε (τ ) exp -i x ±ε/2 λ ∓ (y, τ )dy . (21)
To determine the contraction factor C C Λ,ε to the associated mapping G C2 Λ , the explicit knowledge of the total symbols λ ± is required. However, for a general potential V (x), this is generally impossible. We circumvent this issue by considering an asymptotic expansion {λ ± 1/2-j/2 } +∞ j=0 of the symbols λ ± . To get this estimate, we expand λ ± asymptotically, as the sum of elementary inhomogeneous symbols

λ ± ∼ ±∞ j=0 λ ± 1/2-j/2 ,
and then we truncate this series up to the (p + 1) st term

λ ± ≈ λ ±,p = p j=0 λ ± 1/2-j/2
as proposed in [START_REF] Antoine | Construction, structure and asymptotic approximations of a microdifferential transparent boundary condition for the linear Schrödinger equation[END_REF]. The sign ≈ has to be understood as a truncation of the infinite series by a finite sum, up to a symbol of homogeneity degree -p/2 in τ . In other words, for |τ | large enough, we truncate up to a O(τ -p/2 ). The approximate convergence rate in the elliptic zone is then

C C P,ε ≈ C C Λ,ε ≈ C C,p ε := sup τ ∈Eτ L C,p ε (τ ) (22)
with

L C,p ε (τ ) = exp i ε/2 -ε/2 λ -,p (y, τ ) -λ +,p (y, τ ) dy . ( 23
)
Choosing the principal symbol as λ ± 1/2 = ∓ -τ + V (x), then one gets

λ -,p = -λ +,p , (24)
implying that [START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF] becomes

L C,p ε (τ ) = exp -2i ε/2 -ε/2 λ +,p (y, τ )dy . ( 25
)
Let us remark that [START_REF] Hörmander | Linear Partial Differential Operators[END_REF] does not hold for

λ ± 1/2 = ∓ √ -τ , if V (x) = 0 [6]. A third step
is required to approximate the symbols λ ± 1/2-j/2 , j = 0, ..., p, where 1/|τ | is small (highfrequency regime). For each symbol λ ± 1/2-j/2 , we consider a Taylor expansion up to the order p/2 in 1/|τ |:

λ ±,p ≈ λ ±,p = p j=0 (λ ± 1/2-j/2 ) (1-p)/2 , (26)
where

(λ ± 1/2-j/2 ) (1-p)/2 = λ ± 1/2-j/2 + O(|τ | -p/2 ) .
Notice that Padé's approximants could also be used to approximate λ ± 1/2-j/2 without requiring 1/|τ | small, see [START_REF] Antoine | Absorbing boundary conditions for the onedimensional Schrödinger equation with an exterior repulsive potential[END_REF]. We then define

L C,p ε (τ ) = exp -2i ε/2 -ε/2 λ +,p (y, τ )dy ,
and the associated high-frequency asymptotic convergence rate in the elliptic zone C C,p ε , where

C C P,ε ≈ C C Λ,ε ≈ C C,p ε := sup τ ∈Eτ L C,p ε (τ ). ( 27
)
Let us now set

L ε,1/2-j/2 (τ ) = exp -2i ε/2 -ε/2 λ + 1/2-j/2 (y, τ )dy , L p ε,1/2-j/2 (τ ) = exp -2i ε/2 -ε/2 (λ + 1/2-j/2 ) (1-p)/2 (y, τ )dy . (28)
Therefore, we have

L C,p ε = p j=0 L ε,1/2-j/2 and L C,p ε = p j=0 L p ε,1/2-j/2 . ( 29
)
These preliminary computations lead to the following theorem.

Theorem 2.1. Let V (x) be a smooth potential, and let us assume that the symbols {λ ± 1/2-j/2 } j≥0 are defined as in Proposition 2.2. An asymptotic estimate in the elliptic zone E τ of the contraction factor of the mapping G C2 P defined in [START_REF] Besse | Schwarz waveform relaxation method for one dimensional Schrödinger equation with general potential[END_REF] for the CSWR algorithm (5), is given by

C C P,ε ≈ C C,3 ε = sup τ ∈Eτ L C,3 ε (τ ), (30)
with

L C ε (τ ) ≈ L C,3 ε (τ ) := -τ + V (-ε/2) -τ + V (+ε/2) 1/2 exp -2i ε/2 -ε/2 -τ + V (y)dy (31)
while in the hyperbolic zone H τ the contraction factor is given by the first term in (31).

In addition, one also gets the following approximation of the contraction factor in the elliptic zone, when the symbols are Taylorized for |τ | large

C C P,ε ≈ C C,3 ε = sup τ ∈Eτ L C,3 ε (τ ), (32)
with

L C ε (τ ) ≈ L C,3 ε (τ ) := exp -2ε √ -τ + 1 2 √ -τ ε/2 -ε/2 V (y)dy . (33)
Proof. From (21), we have

êC,± Λ (x, τ ) = h ± ε (τ ) exp -i x ±ε/2 λ ∓ (y, τ )dy .
This implies that

F t G C2 P h + ε , h - ε ≈ exp i ε/2 -ε/2 λ -(y, τ ) -λ + (y, τ ) dy ĥ+ ε , ĥ- ε .
By using Proposition 2.2, one gets

λ ± 1/2 (x, τ ) = ∓ -τ + V (x), λ ± -1 (x, τ ) = ∓ 1 4 V (x) -τ + V (x) , (34)
λ ± 0 = 0 and λ ± -1/2 (x, τ ) = 0. A direct computation leads to L ε,1/2 (τ ) = exp -2i ε/2 -ε/2 -τ + V (y)dy , L ε,-1 (τ ) = exp i 2 ε/2 -ε/2 V (y) -τ + V (y) dy , (35)
with L ε,0 = L ε,-1/2 = 1, and

L ε,-1 (τ ) = -τ + V (-ε/2) -τ + V (+ε/2) 1/2 . ( 36
)
As a consequence, one gets

L C,3 ε (τ ) = -τ + V (-ε/2) -τ + V (+ε/2) 1/2 exp -2i ε/2 -ε/2 -τ + V (y)dy . ( 37
)
We obtain

(λ + 1/2 ) -1 = √ -τ 1 - V (x) 2τ ,
and then

L 3 ε,0 = L 3 ε,-1/2 = L 3 ε,-1 = 1, L 3 ε,1/2 = exp -2ε √ -τ - 1 2 √ -τ ε/2 -ε/2 V (y)dy) (38)
The proof follows from

L C,3 ε = L 3 ε,1/2 .
By neglecting the scattering effects in (6), we have then estimated in the above theorem the approximate rate of convergence for the CSWR method with non-constant potentials.

In the following proposition, we estimate |C C P,ε -C C Λ,ε | assuming that V and V are bounded.

Proposition 2.4. With the same assumptions as in Theorem 2.1, and assuming that V and V are in L ∞ (Ω ± ε ), there exists C(ε, V ) > 0 such that for ε small enough:

C C P,ε ≤ C C Λ,ε + C(ε, V )ε 2 .
This result ensures that if V and V , and/or ε are small enough, the exact contract factor C C P,ε , is close to the approximate one C C Λ,ε computed in Theorem 2.1. This expected result is confirmed numerically in Section 3. Proof. Recall that e C,± P (resp. e C,± Λ ) denotes the solution in L 2 (Ω ± ε ) to (6) (resp. to [START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation[END_REF]). For the sake of notation simplicity, we hereafter omit in the proof the upper index C in e C,± P and e C,± Λ . We now set f ± := e ± Pe ± Λ . Some basic computations show that f ± is solution in

L 2 (Ω ± ε ) to ∂ 2 x -(iτ + V (x)) f ± (x, τ ) = i∂ x λ ± e ± Λ , in Ω ± ε , f ± ± ε/2, τ = 0, at {±ε/2} × R. (39)
As V and V are assumed to be in L ∞ (Ω ± ε ), then e ± Λ (resp. f ± ) for any fixed τ , belongs to

C 1 (Ω ± ε ) (resp. to C 2 (Ω ± ε )).
Multiplying the first equation of (39) by f ± (x, τ ), we get

Ω ± ε f ± (x, τ )∂ 2 x f ± (x, τ )dx = Ω ± ε (iτ + V (x))| f ± (x, τ )| 2 +i∂ x λ ± (x, τ ) e ± Λ (x, τ ) f ± (x, τ )dx ( 40
)
Taking the imaginary part in (40), we have

|τ | Ω ± ε | f ± (x, τ )| 2 dx = Im i Ω ± ε ∂ x λ ± (x, τ ) e ± Λ (x, τ ) f ± (x, τ )dx .
Then for any η > 0, we obtain

|τ | f ± (•, τ ) 2 L 2 (Ω ± ε ) ≤ ∂ x λ ± (•, τ ) L ∞ (Ω ± ε ) η 2 f ± (•, τ ) 2 L 2 (Ω ± ε ) + 1 2η e ± Λ (•, τ ) 2 L 2 (Ω ± ε ) .
By taking η small enough, we deduce that f ± (•, τ ) ∈ L 2 (Ω ± ε ). Then taking the real part in (40) for any η > 0, we get

∂ x f ± (•, τ ) 2 L 2 (Ω ± ε ) ≤ ∂ x λ ± (•, τ ) L ∞ (Ω ± ε ) η 2 f ± (•, τ ) 2 L 2 (Ω ± ε ) + 1 2η e ± Λ (•, τ ) 2 L 2 (Ω ± ε ) + V L ∞ (Ω ± ε ) f ± (•, τ ) 2 L 2 (Ω ± ε ) .
Similarly, multiplying by ∂ 2 x f ± (x, τ ), we have

Ω ± ε |∂ 2 x f ± (x, τ )| 2 dx = -Ω ± ε (iτ + V (x))|∂ x f ± (x, τ )| 2 + V (x) f ± (x, τ )∂ x f ± (x, τ ) dx +i Ω ± ε ∂ x λ ± (x, τ ) e ± Λ (x, τ )∂ 2 x f ± (x, τ )dx ,
and then, for any η > 0,

∂ 2 x f ± (•, τ ) 2 ≤ V L ∞ (Ω ± ε ) + V L ∞ (Ω ± ε ) ∂ x f ± (•, τ ) 2 + V L ∞ (Ω ± ε ) f ± (•, τ ) 2 + ∂ x λ ± (•, τ ) L ∞ (Ω ± ε ) η 2 ∂ 2 x f ± (•, τ ) 2 L 2 (Ω ± ε ) + 1 2η e ± Λ (•, τ ) 2 L 2 (Ω ± ε ) .
By taking η small enough, we deduce that for τ fixed, f ± (•, τ ) ∈ H 2 (Ω ± ε ). Next, we have

F t G C Λ • G C Λ h + ε , h - ε = exp i ε/2 -ε/2 λ -(y, τ ) -λ + (y, τ ) dy h + ε , h - ε as e ± Λ (x, τ ) = h ± ε (τ ) exp -i x ±ε/2
λ ∓ (y, τ )dy .

Now we need to evaluate

F t G C P • G C P h + ε , h - ε
. As f ± is regular in x, by Taylor's expansion and for some ξ ±

x,ε between x and ±ε/2, we have

f ± (x, τ ) = f ± (±ε/2, τ ) + ∂ x f ± (ξ ± x,ε , τ )(x ∓ ε/2) = h ± ε (τ )∂ x g ± (ξ ± x,ε , τ )(x ∓ ε/2)
as f ± (±ε/2, τ ) = 0 and as by definition of f ± , ∂ x f ± is of the form h ± ε (τ )∂ x g ± for some regular function g ± , i.e.

e ± P (x, τ ) = e ± Λ (x, τ ) + h ± ε (τ )∂ x g ± (ξ ± x,ε , τ )(x ± ε/2) .

Now, we can write that

F t G C P • G C P h + ε , h - ε = e - P (ε/2, τ ), e + P (-ε/2, τ ) = h - ε (τ ) exp -i ε/2 -ε/2 λ + (y, τ )dy + ε∂ x f -(ξ - ε/2,ε , τ ), h + ε (τ ) exp -i -ε/2 ε/2 λ -(y, τ )dy -ε∂ x f + (ξ + -ε/2,ε , τ ) = h +,(2) ε (τ ), h -,(2) ε (τ ) .
We iterate once more e ±,( 2)

P (∓ε/2, τ ) = h ±,(2) ε (τ ) exp i ±ε/2 ∓ε/2 λ ∓ (y, τ )dy ∓ ε∂ x f ± (ξ ± ∓ε/2,ε , τ ).
Some additional computations lead to e ±,( 2)

P (∓ε/2, τ ) = h ∓ ε (τ ) exp -2i ε/2 -ε/2 λ + (y, τ )dy ±ε exp -i ε/2 -ε/2 λ + (y, τ )dy ×∂ x f ∓ (ξ ∓ ±ε/2,ε , τ ) -∂ x f ± (ξ ± ∓ε/2,ε , τ ) + O(ε 2 ).
Then as f ± is C 2 , for ε small enough, we obtain e ±,( 2)

P (∓ε/2, τ ) = h ∓ ε (τ ) exp -2i ε/2 -ε/2 λ + (y, τ )dy +ε 2 ∂ 2 x f ∓ (ξ ∓ ±ε/2,ε , τ ) ∓ λ + (0, τ )∂ x f ∓ (ξ ∓ ±ε/2,ε , τ ) + o(ε 2 ) , so that F t G C P • G C P h + ε , h - ε = F t G C Λ • G C Λ h + ε , h - ε + O(ε 2 ), O(ε 2 ) . Let us recall that C C Λ,ε = sup τ exp -2i ε/2
-ε/2 λ + (y, τ )dy . Now using that ∂ x g C,± is bounded for V and V in L ∞ (Ω ± ε ), there exists c(τ, ε, V ) > 0 such that for ε small enough

C C P,ε ≤ sup τ exp -2i ε/2 -ε/2 λ + (y, τ )dy + c(τ, ε, V )ε 2 .
Finally, there exists C(ε, V ) > 0 such that

sup τ F t G C P • G C P h + ε , h - ε = C C P,ε h + ε , h - ε ≤ C C Λ,ε + C(ε, V )ε 2 h + ε , h - ε
This concludes the proof.

Asymptotic estimates of the contraction factor in the q-OSWR algorithm

We consider now the quasi-Optimal Schwarz Waveform Relaxation (q-OSWR) method [START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF] and study its rates of convergence as a function of the order of the transmitting boundary conditions which are used. If we assume that φ ±,(0) ∓ε/2, • and φ ± 0 are some given functions, the q-OSWR algorithm reads as follows for an iteration k ≥ 1:

   P φ ±,(k) = 0, in Ω ± ε × R * + , φ ±,(k) (•, 0) = φ ± 0 , in Ω ± ε , (∂ x + iΛ ±,p )φ ±,(k) ± ε/2, • = (∂ x + iΛ ±,p )φ ∓,(k-1) ± ε/2, • in R * + , (41)
where Λ ±,p = Op(λ ±,p ), for p = 1/2, 0, -1/2... as similarly done in [START_REF] Lorin | Schwarz waveform relaxation domain decomposition methodology for the N-body time-independent and time-dependent Schrödinger equation[END_REF]. Following [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF], we set in Ω ± ε P e O,±,p

P = 0 on Ω ±,p ε × R * + , (∂ x + iΛ ±,p )e O,±,p P ± ε/2, t = h ± ε (t) on {±ε/2} × R * + , (42)
and we introduce the mapping

G O,p P : h + ε , h - ε → (∂ x + iΛ +,p )e O,-,p P ε/2, • , (∂ x + iΛ -,p )e O,+,p P -ε/2, • . (43)
We again want to estimate the asymptotic convergence rate of the q-OSWR domain decomposition method, when considering the system (41) for a transmission operator ∂ x + iΛ ±,p . The q-OSWR method is now associated to the following approximate boundary-value problem on the error function e O,±,p Λ :

(∂ x + iΛ ∓ )e O,±,p Λ = 0, in Ω ± ε × R * + , (∂ x + iΛ ±,p )e O,±,p Λ ± ε/2, t = h ± ε (t), on {±ε/2} × R * + . (44)
Then similarly to [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF], we have Theorem 2.2. Let us assume that V is a smooth one-dimensional space-dependent potential and that λ ± is approximated by P ) defined by (43), for the fixed-point q-OSWR algorithm (41), is given by

λ ± ≈ λ ±,p := p j=0 λ ± 1/2-j/2 ,
C O P,ε ≈ C O,p ε = sup τ ∈Eτ L O,p ε (τ ),
where

L O,p ε (τ ) ≈ c p ε 1 |λ + 1/2 (ε/2, τ )λ + 1/2 (-ε/2, τ )| p+1 L C,p ε (τ ). ( 45
)
In the hyperbolic zone, the frequencies also contribute to the q-OSWR convergence thanks to a contraction factor given by c p ε sup p+1) . In the previous expressions, c p ε is an (ε, p, V )-dependent positive real-valued constant. The principal symbol is given by (9) and L C,p ε (τ ) designates [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]. For a constant potential V , the fixed-point q-OSWR algorithm (41) converges in two iterations for the one-dimensional potential-free case, even without overlap (ε = 0).

τ ∈Hτ |λ + 1/2 (ε/2, τ )λ + 1/2 (-ε/2, τ )| -(
The proof is very similar to the proof of Theorem 2.2 in [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF], where was analyzed the q-OSWR convergence rate for the Schrödinger equation in imaginary time. Except for the symbols, the overall strategy and computational details are identical in both theorems. We here only give the sketch of the proof and again refer to [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF] for more details.

Sketch of Proof.

We first introduce the approximate problem (44) which gives us the approximate representation e O,±,p Λ of the error e O,±,p P . Similarly to the analysis of the CSWR method, the aim is to derive an approximate convergence rate for the q-OSWR method, by introducing an analogous problem to (44). This can be expressed at the symbol level, and we know from (44)

(∂ x + iλ ∓ ) e O,±,p Λ = 0, in Ω ± ε × R * + , (∂ x + iλ ±,p ) e O,±,p Λ ± ε/2, t = h ± ε (t), on {±ε/2} × R * + , (46)
which leads to the following expression

e O,±,p Λ (x, τ ) = α ±,p ε (τ ) exp -i x ±ε/2 λ ∓ (y, τ )dy , (47)
for some functions α ±,p ε . By implementing the transmitting boundary conditions (second equation of system (46)), we get

e O,±,p Λ (x, τ ) = h ± ε (τ ) i(λ ±,p ± ε/2, τ -λ ∓ ± ε/2, τ ) exp -i x ±ε/2
λ ∓ (y, τ )dy .

Now, we consider the mapping

G O,p Λ : h + ε , h - ε → (∂ x + iΛ +,p )e O,-,p Λ ε/2, • , (∂ x + iΛ -,p )e O,+,p Λ -ε/2, • . (48)
By using (47), we get e O,±,p,( 2)

Λ (x, τ ) = λ ∓ ± ε/2, τ + λ ±,p ± ε/2, τ λ ±,p ± ε/2, τ -λ ∓ ± ε/2, τ × h ∓ ε (τ ) i λ ∓,p ∓ ε/2, τ -λ ± ∓ ε/2, τ × exp -i ±ε/2 ∓ε/2 λ ± (y, τ )dy exp -i x ±ε/2 λ ∓ (y, τ )dy .
We recall that λ -= -λ + and λ -,p = -λ +,p , one can show that

F t G O,p Λ • G O,p Λ h + ε , h - ε (τ ) = λ +,p ε/2, τ -λ + ε/2, τ λ + -ε/2, τ -λ +,p -ε/2, τ λ +,p ε/2, τ + λ + ε/2, τ λ + -ε/2, τ + λ +,p -ε/2, τ × exp -2i ε/2 -ε/2 λ + (y, τ )dy h + ε , h - ε (τ ). ( 49
)
From Proposition 2.3, we know that we have the following asymptotic control (for large

| -τ + V | 1) of the remaining λ + (±ε/2, τ ) -λ +,p (±ε/2, τ ) = +∞ j=p+1 λ + 1/2-j/2 (±ε/2, τ ) = 1 (λ + 1/2 (±ε/2, τ )) p +∞ j=p+1 L j+1 =0 F V,-j/2 (λ + 1/2 ) (±ε/2, τ ) = O 1 (λ + 1/2 (±ε/2, τ )) p . (50)
Furthermore, we have

λ + (±ε/2, τ ) + λ +,p (±ε/2, τ ) = 2λ + 1/2 (±ε/2, τ ) + O 1). ( 51
)
Collecting the various estimates, we conclude that

F t G O,p2 Λ h + ε , h - ε (τ ) = O 1 (λ + 1/2 (ε/2, τ )λ + 1/2 (-ε/2, τ )) p+1 × exp -2i ε/2 -ε/2 λ + (y, τ )dy h + ε , h - ε (τ ).
(52) After truncating the symbolic expansion within the exponential term, an approximation rate of convergence for the q-OSWR DDM in both the elliptic zone, is given by

C O,p ε = sup τ ∈Eτ L O,p ε (τ ), (53)
where

L O,p ε (τ ) ≈ c p ε 1 |λ + 1/2 (ε/2, τ )λ + 1/2 (-ε/2, τ )| p+1 L C,p ε (τ ), (54)
where L C,p ε (τ ) is given by ( 25), p ∈ N * , and c p ε is an ε-, p-and V -dependent positive realvalued constant. Unlike the CSWR method, the frequencies from the hyperbolic zone also contributes to the convergence of the q-OSWR method, thanks to the coefficients 1/|τ | p+1 .

x + N + = b + ε/2 and x - j = b -ε/2 + (j -1)∆x, x - N -= a.
The CN scheme is defined as follows: for interior points in (-a, a) and for 1 < j < N ± , and

k ≥ 1 in Ω ± ε φ ±,n+1,(k) j -φ ±,n,(k) j ∆t n = i 2∆x 2 φ ±,n+1,(k) j+1 -2φ ±,n+1,(k) j + φ ±,n+1,(k) j-1 + i 2∆x 2 φ ±,n+1,(k) j -2φ ±,n+1,(k) j + φ ±,n+1,(k) j-1 + i 2 V ± j φ ±,n+1,(k) j + φ ±,n,(k) j , (55)
where

V ± j = V (x ± j ).
The convergence criterion for the Schwarz DDM is given by φ

+,n T ,(k) |Γε -φ -,n T ,(k) |Γε ∞,Γε L 2 (0,T ) ≤ δ Sc , (56)
with δ Sc = 10 -14 ("Sc" for Schwarz) and where n T = T /∆t. When the convergence of the full iterative algorithm is obtained at Schwarz iteration k cvg , one gets the converged global

φ cvg := φ (k cvg ) in Ω a .
Let us remark that (64) is simple to implement since it is a standard CN scheme in each subdomain with homogeneous Dirichlet boundary conditions at the endpoint.

Numerical validation of the CSWR

The CSWR method is based on Dirichlet transmission conditions. At

x + N + = b + ε/2, we impose φ +,n+1,(k) N + +φ +,n,(k) N + = φ -,n+1,(k-1) j 0 +φ -,n,(k-1) j 0
, and j 0 denotes the number of overlapping nodes, i.e. ε = (j 0 -1)∆x. At

x 1 = b-ε/2, we fix φ -,n+1,(k) 1 +φ -,n,(k) 1 = φ +,n+1,(k-1) N + -j 0 +φ +,n,(k-1)
N + -j 0 . Finally at x + 1 = -a and x - N -= a, we set null Dirichlet boundary conditions. The 2 -stability and the second-order accuracy in space and time is straightforwards, as Dirichlet boundary conditions are used. Test 1: potential V = 0. The first test is devoted to two waves propagating in opposite directions from initial time t = 0 to final time T = 1. We follow the experiment given in [START_REF] Antoine | Domain decomposition method and highorder absorbing boundary conditions for the numerical simulation of the time dependent Schrödinger equation with ionization and recombination by intense electric field[END_REF]. The parameters for this first test are the following: a = 10, b = 5/2, with N + = 312 and N -= 193. The size of the overlapping zone is ε = ∆x, corresponding to j 0 = 2. The time step is chosen equal to ∆t = 0.2. The initial data is given by

φ 0 (x) = exp - 1 2 b -a 2 -x 2 exp(ik 0 x) + exp - 3a 8 + x 2 exp(-ik 0 x), with k 0 = 4, see Figure (1) (Left). In Figure (1) (Middle)
, we observe a good agreement between the reconstructed solution using the CSWR algorithm with the reference solution, obtained from a CN scheme on the global domain. From the above analysis (19), we get

C C P,ε = C C Λ,ε = sup L C Λ,ε (τ num ),
where

L C Λ,ε (τ num ) = exp -2iελ + (τ num ) and |τ num | ∈ [1/T, 1/∆t].
The supremum restricted to the elliptic zone, where τ num ∈ [1/T, 1/∆t], exhibits a contraction factor exp -∆x/ √ T < 1, as the supremum in this zone is actually reached at final time T . In Figure [START_REF] Al-Khaleel | Optimized waveform relaxation methods for longitudinal partitioning of transmission lines[END_REF] (Right), we report the numerical convergence rate obtained by the CN scheme, and the theoretical convergence rates (33) in the elliptic zone, and we observe a relatively good agreement. Let us remark that the convergence rate is independent of the time-step ∆t, but depends on the final time T (as the supremum is reached at 1/T ), as reported in Figure [START_REF] Alinhac | Pseudo-differential operators and the Nash-Moser theorem[END_REF], where we observe that the smaller the final time T , the larger the convergence slope. Test 2: potential V = 0. The next test is also devoted to two waves propagating in opposite directions from time t = 0 to T = 1 with a positive space-dependent potential V .

In this case, we expect the following estimate

L C ∆x (τ num ) ≈ τ num -V (-∆x/2) τ num -V (+∆x/2) 1/2 exp -2i∆x ∆x/2 -∆x/2 -τ num + V (y)dy . ( 57
)
The parameters for the second experiment are as follows: a = 10, b = 5/2, with N + = 312, N -= 193 and the time step is ∆t = 0.2. The overlapping region is reduced to two nodes, j 0 = 2, such that ε = ∆x. The initial data is given by

u 0 (x) = exp -4 b -a 2 -x 2 exp(ik 0 x) + exp -4 3a 8 + x 2 exp(-ik 0 x),
with k 0 = 4. The space dependent potential given by V (x) = -20 exp(-5(xb) 2 ). At the interface, we have V (b) ≈ -20. Since V is negative, it is easy to see that the supremum of ( 57) is reached at τ num = 1/∆t, i.e.

C C P,ε ≈ 1/∆t -V (-∆x/2) 1/∆t -V (+∆x/2) 1/2 exp -2i∆x ∆x/2 -∆x/2
1/∆t + V (y)dy .

In Figure (3), we again observe a good agreement between the numerical and theoretical convergence rates. In Figure (4) (Middle) (resp. (Right)), the comparisons of the convergence rates are made between various time steps ∆t (resp. final time T ), showing that the convergence rate, as expected is moderately dependent (resp. independent) of ∆t (resp. T ), as 1/∆t is relatively small compared to V (b) = -20.

In Figure ( 4), we observe that, after one Schwarz iteration with space dependent potential V , the solution is close to convergence unlike the case V = 0.

Numerical validation for q-OSWR

The q-OSWR method involves DtN operators at the subdomain interfaces: with Λ ±,p = Op( λ ±,p ) given by

∂ x φ ±,(k) + i Λ ±,p (x, t, ∂ x , ∂ t )φ ±,(k) = ∂ x φ ∓,(k-1) + i Λ ±,p (x, t, ∂ x , ∂ t )φ ∓,(k-1) , (58)
       Λ +,0 (x, t, ∂ x , ∂ t )φ = -e iπ/4 ∂ 1/2 t φ, Λ +,1 (x, t, ∂ x , ∂ t )φ = -e iπ/4 e iΦ ∂ 1/2 t e -iΦ φ , Λ +,4 (x, t, ∂ x , ∂ t )φ = Λ +,1 φ + 1 4 d n V (x)e iΦ I t e -iΦ φ , (59)
and where the function Φ is defined in the linear case by Φ(x, t) = tV (x). The discretization of the nonlocal time operators are chosen as follows.

∂ 1/2 t f (t n ) ≈ 2 ∆t n k=0 β n-k f k , (60)
I 1/2 t f (t n ) ≈ ∆t 2 n k=0 α n-k f k , (61)
I t f (t n ) ≈ ∆t 2 n k=0 f k , (62)
where the sequence (β n ) n∈N is such that β 0 = 1 and, for n ≥ 0, and (γ n ) n∈N are such that

       (α 0 , α 1 , α 2 , α 3 , α 4 , α 5 , ...) = (1, 1, 1 2 , 1 2 , 3 8 , 3 8
, ...),

β k = (-1) k α k , ∀k ≥ 0, (γ 0 , γ 1 , γ 2 , γ 3 , ...) = (1, 2, 2, ...). (63)
The discretization (60)-( 62) is designed to be consistent with the Crank-Nicolson scheme by using the associated generating function [START_REF] Antoine | Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF]. The semi-discrete q-OSWR-CN scheme for a two-domains decomposition method with a first-order ABC reads as follows

           i φ +,n+1,(k) ∆t + 1 2 ∂ 2 x φ +,n+1,(k) + ∂ 2 x φ +,n,(k) + V (x) 2 φ +,n+1,(k) + φ +,n,(k) = φ +,n,(k) ∆t , in Ω + a,ε , (∂ n ± -e iπ/4 2 ∆t )φ ±,n+1,(k) ±ε/2 = g ∓,n+1,(k-1) ±ε/2 + α ±,n,(k) ±ε/2 -α ∓,n,(k-1) ±ε/2
, φ ±,n+1,(k) = 0, at x = ∓a.

(64) The outwardly directed unit normal vector to Ω ± a,ε is denoted by n ± = ±1. We also set

g ∓,n+1,(k-1) ±ε/2 = ∂ n ± φ ∓,n+1,(k-1) ±ε/2 -e iπ/4 2 ∆t φ ∓,n+1,(k-1) ±ε/2 , α ∓,n,(k) ±ε/2 = e iπ/4 2 ∆t E n ±ε/2 n =0 β n+1-Ē ±ε/2 φ ∓, ,(k) ±ε/2 , E n ±ε/2 = exp -in∆tV ±ε/2 . (65)
Alternatively, a Neumann-to-Dirichlet approximation is possible. For the sake of simplicity, we assume that V = 0, and the scheme reads as follows.

                   i φ ±,n+1,(k) ∆t + 1 2 ∂ 2 x φ ±,n+1,(k) + ∂ 2 x φ ±,n,(k) = φ ±,n,(k) ∆t , in Ω + a,ε , φ ±,n+1,(k) ±ε/2 = -e iπ/4 ∆t 2 n+1 l=0 α l ∂ n + φ ±,n+1-l,(k) ±ε/2 + φ ∓,n+1,(k-1) ±ε/2 +e iπ/4 ∆t 2 n+1 l=0 α l ∂ n + φ ∓,n+1-l,(k-1) ±ε/2
, φ ±,n+1,(k) = 0, at x = ∓a.

(66) We next show a partial stability result for (66). Proposition 3.1. At Schwarz iteration k = 1 and for ε > 0, the semi-discrete scheme (66) is unconditionally L 2 -stable. Moreover, for ε = 0, the semi-scheme (66) is unconditionally L 2 -stable at any Schwarz iteration.

Proof. We multiply i) -i( φ+,n+1,(k) + φ+,n,(k) )/2 (resp. -i( φ-,n+1,(k) + φ-,n,(k) )/2), ii) integrate on Ω + ε (resp. on Ω - ε ) and take the real part:

1 2∆t φ +,n+1,(k) 2 L 2 (Ω + ε ) = 1 2∆t φ +,n,(k) 2 L 2 (Ω + ε ) + Re i φ+,n+1,(k) + φ+,n,(k) 2 ∂ n + φ +,n+1,(k) + φ +,n,(k) 2 +ε/2 Set φ +,n+1/2,(k) +ε/2 = φ +,n+1,(k) +ε/2 + φ +,n,(k) +ε/2 2 , φ -,n+1/2,(k) -ε/2 = φ -,n+1,(k) -ε/2 + φ -,n,(k) -ε/2 2 so that 1 2∆t φ +,n+1,(k) 2 L 2 (Ω + ε ) = 1 2∆t φ +,n,(k) 2 L 2 (Ω + ε ) + Re i φ+,n+1/2,(k) ∂ n + φ +,n+1/2,(k) +ε/2 and 1 2∆t φ +,n+1,(k) 2 L 2 (Ω + ε ) = 1 2∆t φ +,0,(k) 2 L 2 (Ω + ε ) + n l=0
Re i φl+1/2,(k) ∂ n + φ l+1/2,(k) +ε/2

Similarly, we have

1 2∆t φ -,n+1,(k) 2 L 2 (Ω - ε ) = 1 2∆t
φ -,0,(k) 2 In other words for all k ≥ 1, we get Re i φ±,l+1/2,(k) ∂ n ± φ ±,l+1/2,(k) 0 ≤ 0 . This ensures the L 2 -stability of the semi-discrete scheme. The proof of stability at any Schwarz iteration when > 0 requires additional assumptions and is not presented in this paper.

Based on this discretization, we now validate the theoretical rate of convergence for the q-OSWR method. Let us define T = T 0 /2 , with = -2, -1, 0, 1, 2. According to the above analysis, we have C In Figure 5, we compare the theoretical and numerical slopes for p = 0. We observe that the numerical convergence slope is close to the theoretical one, i.e.log(2)/2. To validate the theoretical slopes for larger values of p, it is necessary to develop a much heavier numerical machinery not presented in this paper.

Conclusion

The paper was devoted to the analysis of the rate of convergence of the Schwarz Waveform Relaxation domain decomposition methods applied to the linear Schrödinger equation with space dependent potential. Approximations of the theoretical convergence rate have been stated for both the CSWR and the q-OSWR algorithms by using pseudo-differential. In addition, they were also numerically validated. In a forthcoming paper, we will analyze the convergence rate of SWR in higher dimensions and for more subdomains.

Proposition 2 . 3 .

 23 Let us define the class of symbols S

  with p ∈ N. An asymptotic estimate in elliptic region of the contraction factor C O P,ε of the mapping G O,p2 P (with G O,p2 P := G O,p P • G O,p

Figure 1 :

 1 Figure 1: (Left) Initial data amplitude. (Middle) CSWR: reconstructed converged solution after iteration k until convergence at final time T . (Right) Theoretical and numerical convergence rate comparison.

Figure 2 :

 2 Figure 2: Convergence rates for different final times T .

Figure 3 :

 3 Figure 3: (Left): Reconstructed solution, and reference solution with space dependent potential V . (Right) Theoretical and numerical convergence rate comparison.

2 Figure 4 :

 24 Figure 4: (Left) Reconstructed solution after 1 Schwarz iteration, and reference solution at final time T = 1. (Middle) Convergence rates with different time-steps ∆t and space-dependent potential V for CSWR. (Right) Convergence rates with different final times T , and space-dependent potential V for CSWR.

L 2 2 and ζ = e iπ/4 ∆t 2 .±+

 222 Re i φ-,l+1/2,(k) ∂ n -φ -,l+1/2,Re i φ±,l+1/2,(k)∂ n ± φ ±,l+1/2,(k) ±ε/The transmission conditions can now also read φ Now at k = 0, we assumed that φ using Theorem 4 from[START_REF] Antoine | Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF], we deduce Re i φ±,l+1/2,(l) ∂ n ± φ ±,l+1/2,schemes in Ω + ε and Ω - ε , which are then purely non-reflecting, are then unconditionally L 2 -stable if k = 1.Note that for ε = 0, we directly get from the transmission conditions that φ ±,n+1,(k) 0+ ζ n+1 l=0 α l ∂ n ± φ ±,n+1-l,(k) 0 = φ ∓,n+1,(k-1) 0 + ζ n+1 l=0 α l ∂ n ± φ∓,n+1-l,(k-1) 0 .

  /2 , where p ∈ N designates the order of approximation of the "transparent" operator Λ + and c p ε is an unknown (ε, p)-dependent constant. In order to numerically validate this result, we first observe that logC (OSWR) p, ≈ log C (CSWR) + log c p ε T (p+1)/2 0 -(p + 1) log(2)/2.We next numerically evaluate| log C (OSWR) p, log C (CSWR) |as a function of . For p = 0, the transparent transmission operator is approximated by ∂ 1/2 t the expected slope islog(2)/2, when plotting log C (OSWR) p, as a function of . The numerical data are as follows: N + = 641, N -= 385, a = -20, ∆t = 0.1, T 0 = 1, ε = ∆x, b = 4 and k 0 = 4, and φ 0 (x) = exp -4 bik 0 x), and V (x) = 5 exp -5(xb) 2 .

Figure 5 :

 5 Figure5: Theoretical and numerical convergence rate for q-OSWR with p = 0.

Acknowledgments.

X. Antoine thanks the support of the french ANR grants "Bond" (ANR-13-BS01-0009-01) and ANR-12-MONU-0007-02 BECASIM (Modèles Numériques call). E. Lorin thanks NSERC for the financial support via the Discovery Grant program.

Convergence of the CSWR and q-OSWR algorithms

We recall here two convergence theorems for the CSWR (5) and q-OSWR (41) algorithms when V is constant, and which were proven in [START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF]. Unlike, Theorem 1.1, their proof relies on the calculation of the contraction factor, as proposed above.

Theorem 2.3. Let the initial guess (h -, h + ) be in (H 3/4 (0, T )) 2 and such that h -(0) = φ 0 (-ε/2) and h + (0) = φ 0 (ε/2). Let V be a real constant and suppose that τ max > 0 is such that F(e -t h ± ) vanishes outside of [-τ max , +∞) with g + = h +φ + (-ε/2, •) and g -= h -φ -(ε/2, •). Then the iterates (φ -,(k) , φ +,(k) ) of the algorithm [START_REF] Antoine | Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF]

to the solution of (1).

The above theorem describes the convergence of the overlapping Schwarz waveform relaxation algorithm, stating that the convergence is at least linear and that it depends on the size of the overlapping region. For the CSWR, if the potential V is constant, the approximate and the exact contraction factors coincide, that is

The q-OSWR convergence theorem reads as follows:

Theorem 2.4. Let V be a real constant. The algorithm (41) converges to the solution φ of (1) in two iterations for any initial guess φ +,(0) and φ -,(0) , independently of the size of the overlap (-ε/2, ε/2) for ε > 0 if and only if

When V is constant, the OSWR algorithm then converges in two iterations, which corresponds to a contraction factor equal to 0. When V is non-constant, we expect a similar result. Thanks to the presented analysis (at least at high frequencies), the contraction factor can be shown to be arbitrarily small in the elliptic zone. However, notice that a rigorous proof for V non-constant is yet to be established (work in progress).

Numerical validation

Numerical approximation of the interior scheme

In the one-dimensional case, we consider two bounded subdomains

Finally, b ± ε/2 denotes the position of the interfaces. Notice that in the above analysis, we have taken b = 0. The numerical scheme that we use is a second-order unconditionally 2 -stable Crank-Nicolson (CN) scheme [START_REF] Antoine | Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF][START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]. We define by ∆x the space step, and by ∆t n the time-step at iteration n. We denote by φ ±,n,(k) j the solution in Ω ± ε of the scheme at Schwarz iteration k, at time iteration n and at x ± j , where the nodes {x ± j } 1≤j≤N ± are defined by: x + j = -a + (j -1)∆x,


Saif Al-Kuwari

Stephen Wolthusen

FORENSIC TRACKING AND MOBILITY PREDICTION IN VEHICULAR NETWORKS

Keywords: Vehicular networks, passive tracking, mobility prediction

Vehicular networks have attracted significant attention in the context of traffic control and management applications. However, vehicular networks also have important applications in crime prevention and criminal investigations. This paper presents a system for passively tracking a target vehicle whose driver is assumed to be a "person of interest." The tracking system relies on the dynamic recruitment of neighboring vehicles of the target as agents. A mobility prediction algorithm is used to probabilistically predict the target's future movement and to adjust the tracking process. Combining agent-based tracking and mobility prediction enables a target vehicle to be passively localized and tracked in an efficient manner.

Introduction

Due to the heightened interest in driver safety and infotainment applications, vehicular networks along with other mature wireless technologies will be widely deployed in future automobiles. This has motivated the research community to investigate various aspects of vehicular networks.

The main motive behind the emergence of vehicular networks is driver safety. However, these networks have important applications in crime prevention and criminal investigations in which law enforcement agencies must track the movements of "persons of interest."

Vehicular network algorithms differ from those encountered in conventional mobile networks. Vehicular network nodes move more rapidly than the nodes in other mobile networks. However, vehicular nodes are characterized by somewhat limited freedom of motion since their move-ments are constrained by roadways and traffic regulations. This introduces additional challenges for vehicular tracking algorithms because they must efficiently adapt to the peculiarities of vehicular movements; but this also improves their mobility predictions because nodes move along predefined paths.

This paper presents a passive vehicular tracking system that relies on the dynamic recruitment of neighboring vehicles as agents. A mobility prediction algorithm is used to probabilistically predict the target's future movement and to adjust the tracking process. Combining agentbased tracking and mobility prediction enables a target vehicle to be localized and tracked both efficiently and clandestinely.

Related Work

Boukerche, et al. [START_REF] Boukerche | Vehicular ad hoc networks: A new challenge for localization-based systems[END_REF] discuss localization techniques that can be used in vehicular networks along with their practical implications. The authors demonstrate that most localization techniques suffer from inherent inaccuracies that may not be acceptable for vehicular-based applications that require precise location information. In such situations, the best solution is to use data fusion where the results of multiple localization techniques are combined to increase accuracy.

While localization techniques for vehicular networks are usually GPSbased, not all vehicles are equipped with GPS devices. Also, GPS-based techniques are useless when GPS signals are not available (e.g., inside tunnels). Benslimane [START_REF] Benslimane | Localization in vehicular ad hoc networks[END_REF] addresses this problem in an extension to the ODAM messaging dissemination protocol, enabling vehicles that are not GPS-capable to be localized.

A popular tracking technique that is well suited to vehicular scenarios involves map-matching tracking, which attempts to match a node's actual location (raw) data to maps. Barakatsoulas, et al. [START_REF] Brakatsoulas | On map-matching vehicle tracking data[END_REF] present several such algorithms that exploit vehicular trajectory information. However, this approach is not suitable for applications that require realtime tracking. Other tracking techniques such as installing tracking tags on a target vehicle [START_REF] Mansell | Vehicle Tracking and Security System[END_REF] also exist, but these solutions are not ad hoc and require additional preparation and overhead.

Vehicular Networks

Generally, vehicular networks are based on ad hoc infrastructures and are, therefore, referred to as vehicular ad hoc networks (VANET), a subclass of mobile ad hoc networks (MANETs). The U.S. Federal Communications Commission (FCC) has allocated a 75 MHz spectrum in the 5.9 GHz band (5.850 GHz to 5.925 GHz) for vehicular communications, which can either be vehicle to vehicle (V2V) or vehicle to infrastructure (V2I) communications. V2V depends on an ad hoc infrastructure where vehicles directly exchange information like accident and congestion warnings. V2I, on the other hand, assumes the presence of pre-installed road components that vehicles communicate with in order to retrieve information. However, the cost of installing these components often limits V2I applications.

Communication in vehicular networks is based on the conventional IEEE 802.11 wireless standard (WiFi) or 3G via CDMA (code division multiple access) technology. WiFi is simpler and less expensive to deploy, but has lower reliability because it was not designed to operate in environments with rapid movements. In contrast, 3G is more relevant to vehicular communications, but is more expensive and difficult to deploy because of its centralized architecture [START_REF] Luo | A Survey of Inter-Vehicle Communication[END_REF]. In this paper, we assume that vehicles communicate using WiFi, which is a more widely deployed solution. Note that WiFi has some efficiency limitations, especially because it operates in the half-duplex mode (i.e., two nodes cannot communicate with each other simultaneously).

Mobility in Vehicular Networks

Mobility models represent the motion behavior of vehicles. Mobility models designed to simulate the movement of nodes in MANETs can be used to simulate vehicular movements in VANETs and most vehicular simulation approaches employ these models. However, MANET models have proved unreliable and unrealistic, which has motivated the design and development of mobility models targeted for vehicular networks.

Vehicular mobility models are categorized as microscopic, mesoscopic and macroscopic models. Microscopic models describe the dynamics of the individual movement of each vehicle and its interaction with other vehicles (this may require excessive computational resources for largescale simulations). Macroscopic models simulate the characteristics of the roadways (motion constraints) and the flow of vehicles, but do not consider individual vehicle movement, which impacts their degree of realism. Interested readers are referred to [START_REF] Baumann | Towards realistic mobility models for vehicular ad hoc networks[END_REF] for a discussion of macrosopic models versus microscopic models. Mesoscopic models belong to a new "hybrid" class of vehicular mobility models that combine features of microscopic and macroscopic models. In this paper, we adopt the microscopic intelligent driver motion (IDM) model [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF] along with two extensions, IDM-IM (IDM with intersection management) and IDM-LC (IDM with lane changes).

Vehicle Localization

Conventional localization techniques used in cellular and sensor networks can also be applied to vehicular networks, albeit with slightly degraded accuracy. These techniques involve measuring the received signal strength (RSS), time of arrival (TOA) or time difference of arrival (TDOA). Since we adopt a passive tracking approach, only RSS is relevant because, unlike TOA and TDOA, it does not require a fully synchronized network and active communication with the target.

The RSS technique involves measuring the strength of the signals emitted from the target vehicle in order to estimate the distance to the target. RSS measurements are nonlinear due to radio propagation effects; nevertheless, our algorithm adopts a free space radio propagation model, which assumes a direct line of sight path between vehicles. This assumption does not severely degrade the accuracy of the localization process because it is most likely that the only obstacles would be the moving vehicles and their effects can be eliminated by averaging the measurements. Note that vehicular acceleration also affects RSS measurements, but we assume that this effect can be mitigated by adding noise in the computations.

The localization of a target vehicle C s involves three of its neighboring vehicles. However, only two of these neighbors (tracking agents) track the target; the other (localizing agent) is retired immediately after the localization is completed. We denote the tracking agents as C a (main) and C b (backup) and the localizing agent as C l . We assume that the three agents can localize themselves using GPS, which is fairly common in current vehicles.

The localization algorithm has four steps:

Step 1 (Recruitment): The tracking agents are recruited by randomly requesting RSS readings from nodes around the target.

The nodes with the highest readings are recruited. Details about the initial and subsequent recruitment process are provided in Section 7.

Step 2 (Distance Measurement): The distances between the recruited agents and the target are estimated using the Friis equation [START_REF] Friis | A note on a simple transmission formula[END_REF]:

d = P t A r /4πP r
where P r and P t are the received and transmitted power, respectively; and A r is the effective area of the receiving antenna.

Step 3 (Trilateration): The location of the target is estimated using a geometric transformation to find the intersection point of the three circles centered at the agents that have radii equal to the respective distances of the agents from the target as measured in

Step 2.

Step 4 (Retirement): The localizing agent C l is retired.

The localization algorithm can also be used in conjunction with the prediction algorithm (Section 6) to determine when the next recruitment is required. This assumes that the agent can measure its own speed, which is possible because it is equipped with a GPS device. The average speed is estimated by dividing a particular distance by its corresponding traversal time.

Mobility Prediction

Mobility prediction is an important aspect to consider when tracking vehicles. Generally, mobility prediction in vehicular networks is slightly simpler in VANETs than MANETs because vehicular mobility is restricted to roadways and the pace of movement is usually limited to the maximum allowable speed of the roadway. The motion of mobile nodes in MANETs, on the other hand, is often unconstrained, although it may be restricted by terrain characteristics and physical limitations. Note also that vehicular movements have higher accelerations.

This section presents a vehicular mobility prediction algorithm that probabilistically predicts the near-future movement of the target vehicle based on its current location and estimated speed, assuming that the target has already been localized. The algorithm incorporates: (i) time prediction, which estimates the time elapsed before the target reaches the next intersection; and (ii) direction prediction, which identifies the direction that the target will most likely take after passing the intersection.

Time Prediction

Predicting the time taken for a vehicle to reach the next intersection requires an estimate of the speed of the vehicle. The estimation procedure is illustrated in Figure 1. Vehicle C a estimates the speed of C s by making two RSS measurements at times t 1 and t 2 . For reasons of simplicity, it is assumed that the distance between the horizontally-aligned vehicles in a roadway is I, which is easily obtained. However, it is very unlikely that vehicles C a and C s will be aligned perfectly. Consequently, when the first RSS measurement (RSS 1 ) is taken, the vertical distance ahead or behind C a is calculated to be perfectly adjacent to C s . The position of C a is then adjusted so that it is aligned with C s using:

N = (d 1 (t 1 )) 2 -I 2
as shown in Figure 1(a). Thus, the position C a becomes C a+N as shown in Figure 1(b).

After the position of C a is adjusted, the second RSS measurement (RSS 2 ) is taken at time t 2 . However, the distance obtained from RSS 2 must be adjusted by recalculating it as if it was measured by C a+N , which accounts for the extra distance N .

Figure 1(c) illustrates this process. First, the distance M is calculated as:

M = (d 2 (t 2 )) 2 -I 2 -N and thus, d 2+N (t 2 ) = M 2 + I 2 .
Next, the speed of C s is computed as:

v s = F + E t 2 -t 1 (1)
where

F = (d 2+N t 2 )) 2 -I 2 and E = v a (t 2 -t 1 ) assuming that the speed v s of C a is known (v s is estimated easily).
After the speed is estimated, the time taken for C s to reach the next intersection is calculated using the traditional distance equation:

t z = n -d x + ϵ/v s
where d x is the distance the vehicle has covered on a roadway of length n and can be calculated relative to the location of C a ; ϵ is a uniformly distributed random variable introduced to compensate for the acceleration/deceleration error margin; and v s is the estimated speed of vehicle C s according to Equation 1.

It is very unlikely that vehicles C a and C s have perfect horizontal alignment; thus, we assume that one vehicle is leading or lagging the other. In Figure 1, vehicle C s is leading C a and is moving faster than C a , so it will continue to lead vehicle C a . Note that the algorithm accommodates other scenarios where C a is the leading vehicle. This is because the absolute velocity of C s is measured without considering its relation to the location and speed of C a .

Since localization is performed before the prediction algorithm is executed, the lane that the target occupies is known. The discussion above assumes that C a and C s are not in the same lane, but the probability of them being in the same lane is equally likely. If this is the case and if C a makes the first RSS measurement (RSS 1 ) and later makes the second RSS measurement (RSS 2 ), then the distance C s has moved during the period between these two measurements is simply the difference between RSS 2 and RSS 1 . Moreover, if the vehicles are not in the same lane, it does not matter which vehicle is in which lane. In Figure 1, it happens that C a is in the left lane and C s is in the right lane, but the algorithm still works if the lane positions are reversed. However, to simplify the design of the algorithm, we assume that the leading entity will continue to lead during the period between the RSS 1 and RSS 2 measurements. This is a reasonable assumption because the time period is usually short enough to preserve the leading status.

Direction Prediction

Predicting the direction that a target vehicle will take through the next intersection is more complex and bears a probabilistic distribution. Also, it presupposes that the target vehicle has been accurately localized in order to identify the lane it occupies. We assume that all drivers adhere to the following basic traffic rules as illustrated in Figure 2:

Rule 1: For a vehicle to turn right at an intersection, it must be in the right lane of the roadway leading to the intersection.

Rule 2: For a vehicle to turn left at an intersection, it must be in the left lane of the roadway leading to the intersection. Rule 3: For a vehicle to go straight ahead at an intersection, it can be in the right lane or in the left lane of the roadway leading to the intersection. Rule 4: A vehicle must position itself in the appropriate lane based on its intended direction and Rules 1-3 around 300 meters before the next intersection or halfway through the roadway leading to the next intersection, whichever comes later. A lane change beyond this point for the purpose of accelerating the pace (e.g., to overtake a vehicle) is not permitted.

It is reasonable to assume that the target vehicle will obey these traffic rules because the target (who is either a suspect or criminal) will most likely not want to attract attention. Obeying Rules 1-4 forces a vehicle to start the lane changing process when reaching approximately n/4 of a roadway of length n or 400 meters away from the intersection, whichever is later. This should allow sufficient time for the vehicle to complete the lane changing process without violating Rule 4. We call the region between n/4 and n/2 (or similarly, the 100 meters between 400 meters to 300 meters away from the intersection) a "critical area" and apply the following predictions:

IF C L i CriticalArea --------→ C R i T HEN P r(C i → R) = 80% IF C R i CriticalArea --------→ C L i T HEN P r(C i → L) = 70%
Based on these predictions, if a particular vehicle C i occupied the left lane (C L i ) and shifted to the right lane as it passed through the critical area, the algorithm infers that the shift was imposed by traffic regulations. This means that the vehicle is very likely going to turn right (R) at the next intersection because: (i) the probability of turning left (L) from a right lane (where the vehicle has shifted) is 0; and (ii) if the vehicle intended to go straight ahead (A), it would most likely have stayed in its original lane without taking the overhead of a lane change as it is permitted to go straight ahead from either lane.

A similar argument applies to a vehicle that shifted from the right lane to the left lane as it passed through the critical area. The difference is that it is possible that the vehicle shifted from the right lane to the left lane to increase its pace of the movement while actually intending to go straight ahead. This decreases the probability of turning left for this particular course of action.

On the other hand, if the vehicle did not change its lane as it went through the critical area, it has a probability of 50% of taking either of the two permitted directions depending on its current lane. Also, it has a 0% probability of taking the banned direction (i.e., left for the right lane and right for the left lane).

An interesting extension to this direction prediction technique (which is not currently implemented in the algorithm) is to integrate points of interest [START_REF] Engelhart | A spatial analysis of mobility models: Application to wireless ad hoc network simulation[END_REF] corresponding to frequently-visited locations such as grocery stores, banks, restaurants and offices. Such locations can significantly influence the probability distribution of vehicular mobility. For example, we know that the probability of a vehicle in the right lane turning right is (currently) the same as the probability of the vehicle going straight ahead if it had not shifted to the right lane as it passed through the critical area. However, that fact that a point of interest is located on the right could increase the probability of the vehicle turning right instead of going straight ahead. Note that the specific time of day can significantly influence the effect of a point of interest on the prediction probabilities (e.g., an office building will most likely be a popular point of interest only during the daytime).

Tracking

In forensic and law enforcement applications, vehicle tracking is usually required to be passive so that the driver of the target vehicle is not aware of the tracking process. This passivity requirement potentially eliminates more active tracking systems such as those involving vehicle telematics.

In our tracking scenario, a group of trackers C t(i) (i = 1, 2, ..., n) consisting of n privately-connected vehicles (e.g., police patrols) recruit a main tracking agent C a and a backup tracking agent C b from the public to track a target vehicle C s . We assume that the recruiting C t(i) is initially located in range of C s , which enables it to recruit suitable C a and C b possibly by visual estimation. During this initial stage, the C t(i) also recruits an additional localization agent C l which, along with the tracking agents, localize the target as discussed in Section 5. The recruiting C t(i) supplies an address list of the other valid C t(i) to all the recruited agents.

The localizing agent is retired after the target is localized, leaving C a to be responsible for the rest of the tracking. C a is backed up by C b , whose responsibility is to take over the tracking process if C a suddenly fails. C a is also responsible for sending location updates of C s to a C t(i) whenever localization is triggered. C a must ensure that its backup agent C b is within range of C a and the target at all times by regularly monitoring its RSS and probing it for the target's RSS; otherwise it has to recruit another C b .

When tracking is initiated, the C t(i) creates a tracking table T rackT B, which stores records of the target's movements as received from C a . This table is synchronized with all other C t(i) as soon as an update is received by a C t(i) . When a localization update is available, C a searches its range for a C t(i) to update it; if multiple C t(i) are found, C a randomly chooses one of them; this potentially minimizes a traffic analysis attack by a third observer who observes the traffic between the agents and the trackers. If no C t(i) is available, C a creates a temporary tracking table T empT B and accumulates information while regularly probing for a valid C t(i) . When a C t(i) is found, C a transfers its T empT B to the C ti which, in turn, merges it with its copy of T rackT B and synchronizes it with the other C t(i) .

After running the prediction algorithm, C a estimates how long it will be able to track the target (called the "alive period") and schedules the next localization for near the expiration of this period. During the alive period, C a keeps observing the RSS measurements from C s . If a specific lower threshold of RSS is reached or the predicted alive time is near expi-Algorithm 1 Vehicular Tracking Algorithm end if 23: until Tracking Expires ration, C a initiates a "probe process," which involves sending requests to the neighboring vehicles, supplying the C s address and asking for their RSS measurements if they can receive transmissions from C s . Note that the vehicle with the strongest RSS is recruited as the localization agent C l . Furthermore, a C t(i) may at any time request C a to localize C s , in which case C a executes the probe process to recruit C l and works with C b to localize C s .

If C a had to maintain a T empT B, but attempted to recruit another C a (due to the expiry of the alive period or weak RSS observations from C s ) and there is still no C t(i) in range, the old C a recruits a new C a and hands off the tracking process to it along with the T empT B. The T empT B is maintained by the new C a until a valid C t(i) is found upon which time the old C a retires itself. Details of the tracking process are presented in Algorithm 1.

Simulation Results

A realistic vehicular environment was created using the VanetMo-biSim simulator [START_REF] Harri | VanetMobiSim: Generating realistic mobility patterns for VANETs[END_REF]; vehicular mobility traces were generated based on IDM-LC model. The ns-2 simulator was used to run the tracking scenarios.

The prediction algorithm was evaluated by executing it over a number of different scenarios and estimating the time taken for vehicles to reach the next intersection. These times were compared with the actual times taken to reach the intersection as reported by the network simulator. The simulation scenarios adopted a set of interconnected roadways with different lengths and speed limits. Additionally, mobility traces with different node densities were generated in order to analyze the effect of node density on the accuracy of the prediction algorithm.

Figures 3 and4 present the results of the prediction algorithm for node densities of 10 and 100 nodes with the simulations running for 1,000 seconds. Clearly, there is a variation in performance. The main factors influencing the precision of the prediction algorithm are: Node Density: As seen in Figures 3 and4, the accuracy of the algorithm is moderately affected by increased node density. This is expected in a real-world environment -the more vehicles present on a roadway, the harder it is to accurately predict their movement due to their irregular acceleration/deceleration behaviors.

Roadway Length: The length of the roadway affects the prediction algorithm because longer roadways allow for more fluctuations in acceleration, which can be difficult to model. Highways have very different characteristics than city roadways in terms of traffic flow.

Speed Limit: Roadways with higher speed limits also allow for increased acceleration fluctuations, which affect the accuracy of the algorithm.

Duration of the Prediction: Other than the macroscopic factors above, the duration that the prediction process is applied to the target is an important factor. The prediction algorithm measures the time taken by the target to pass a specific distance and then adds noise to the measurement to compensate for future acceleration/deceleration before calculating the predicted time. However, the accuracy of the algorithm is improved when the target is observed for a longer interval. This is illustrated in Figures 3 and4, where a longer prediction interval of 20 seconds (i.e., the target's movement is observed for 20 seconds) yields a better prediction than a shorter interval. Another way to improve prediction accu-racy is to take multiple observations at different times and compute the average.

Location of the Prediction: Despite the fact that the prediction process can take place at any point on the roadway, the point where the prediction process takes place can be important, especially for highways that typically have higher speed limits than city roadways. For example, it is not advisable to execute the prediction process at the beginning of a highway. Since vehicles just start accelerating at this point, their accelerations are not representative of the rest of the journey.

Based on the above factors, it is clear that the tracking environment has an unavoidable influence on the tracking process. Carefully modeling the macroscopic features of the environment is, therefore, critical to improving tracking accuracy. Information about the general structure and layout of roadways can be obtained from online databases such as TIGER [START_REF] Census | Topologically Integrated Geographic Encoding and Referencing[END_REF]. Other important road characteristics (e.g., density) can be modeled and estimated based on empirical observations.

Conclusions

Vehicular tracking has important applications in crime prevention and criminal investigations. The algorithms presented in this paper implement the passive localization, tracking and prediction of the movement of a target vehicle. While there has been limited research related to passive vehicular tracking, we argue that the passivity requirement is very important for forensic and law enforcement purposes. Note, however, that the agent-based approach may face barriers in some jurisdictions where randomly recruiting tracking agents may be restricted by law.

Finally, we note that the RF measurements used in the algorithms include an unavoidable error margin. However, estimating the measurement parameters can provide an acceptable error margin. This is because in vehicular tracking it is only important to detect the presence of the target and partially track it to estimate its most likely trace in the corresponding roadway. Inferences can then be made if the target halted en route (e.g., to commit a crime) by observing the estimated time delay along each roadway.

Figure 2 .

 2 Figure 2. Traffic rules at an intersection.

1 :repeat 6 : 7 : 8 : 9 :

 16789 C t(i) ⇐ Trackers {identify n trackers} 2: C s ⇐ Target {identify the target} 3: Recruit(C a , C b ) {recruit tracking agents} 4: Recruit(C l ) {recruit localization agent} 5: Localize(C s ) {C a , C b , C l Localize C s } Retire(C l ) Alive ← P redict(C s ) {find alive period based on the prediction algorithm} while currentT ime ≤ Alive -ω do 10: if C a .RSS(C s ) ≤ threshold then 11: break {if RSS from C s ≤ threshold, break} 12: a .Probe {search for agents} 15: if C a .RSS(C s ) ≤ threshold then 16: if C b .RSS(C s ) ≤ threshold then 17: Recruit(C a , C b , ) {recruit new C a and C b } a ) {recruit new C a } 20:

  Recruit(C l ) {recruit C l } 22:

Figure 3 .

 3 Figure 3. Simulation results (10 node scenario).

Figure 4 .

 4 Figure 4. Simulation results (100 node scenario).