Stationary Schrödinger equation in the semi-classical limit: numerical coupling of oscillatory and evanescent regions

This paper is concerned with a 1D Schrödinger scattering problem involving both oscillatory and evanescent regimes, separated by jump discontinuities in the potential function, to avoid “turning points”. We derive a non-overlapping domain decomposition method to split the original problem into sub-problems on these regions, both for the continuous and afterwards for the discrete problem. Further, a hybrid WKB-based numerical method is designed for its efficient and accurate solution in the semi-classical limit: a WKB-marching method for the oscillatory regions and a FEM with WKB-basis functions in the evanescent regions. We provide a complete error analysis of this hybrid method and illustrate our convergence results by numerical tests.


Introduction
This paper deals with the design, error analysis, and numerical study of an asymptotically correct scheme for the numerical solution of the stationary Schrödinger equation in one dimensional scattering situations: where 0 < ε 1 is a very small parameter and a(x) a piecewise (sufficiently) smooth, real function. On the one hand, for a(x) > 0, the solution is highly oscillatory, with the small (local) wave length λ(x) = 2πε √ a(x) . On the other hand, for a(x) < 0, the wave function ψ is (essentially) exponentially de/increasing, typically exhibiting a thin boundary layer with thickness of the order O ε √ |a(x)| . A key aspect of this paper is that a(x) takes both signs. Hence, we have to cope with a classical multi-scale problem, combining different types of arduousnesses and multi-scale behaviours. Numerically, we aim at recovering these fine structures of the solution, however without using a fine spatial grid. To this end we shall develop a (non-overlapping) domain decomposition method (DDM) to separate the oscillatory and evanescent regions, as they require very different numerical approaches. This DDM allows to recover at the continuous level the exact analytical solution in a single sweep (against the direction of the incoming plane wave) with appropriate interface conditions and a final scaling.
The study of such multi-scale problems is very challenging from a theoretical as well as numerical point of view. In both situations (or regions) a classical discretization (like in [10,11]) requires a very fine mesh in order to accurately resolve the oscillations and boundary layers. However, with a step size requirement of h = O(ε), standard numerical methods would be very costly and inefficient here.
Concerning the oscillatory case, several WKB-based numerical schemes (named after the physicists Wentzel, Kramers, and Brillouin) have been developed and analysed in the last decade. Their goal is to use a coarse spatial grid with step size h > λ (see Figure 1), reducing the limitation to at least h = O( √ ε). For marching schemes we refer to [2,13,17], whereas a finite element method (FEM) using oscillatory WKBbasis functions was introduced in Sect. 2 of [19] and in [20]. This FEM approach has the disadvantage that it requires a non-resonance condition between the mesh-size h and the wave length λ of the solution. By contrast, this restriction is not necessary in the above mentioned marching schemes. Numerical approaches for the evanescent regime (as ε → 0) have been considered much less, so far. We refer to Sect. 3 in [19] for the formulation of a WKB-based multiscale FEM-scheme, including its numerical coupling to the oscillatory region (also based on a FEM). A detailed comparison (w.r.t. accuracy and simulation time) of this WKB-based FEM with a standard FEM was carried out in Sect. 3.1 of [5]. But a numerical analysis has, to our knowledge, not been carried out yet. This paper also aims at closing this gap. In this evanescent regime the problem (1.1) is elliptic, and for the example of a = const., a solution is given by a linear combination of the basis-functions exp(± √ |a| ε x). Hence this region must be considered as a boundary value problem (BVP) and solved e.g. by a finite difference or a FEM method. A In standard numerical methods highly oscillating solutions require a very fine mesh to capture the oscillations. However, with the analytic pre-processing of our WKB-marching method an accurate solution can be obtained on a coarse grid (dots). Plotted is the solution Re[ψ(x)] of (1.1) with ε = 0.01, h = 0.125, and a = (x + 1 2 ) 2 reformulation as an initial value problem (IVP) and the use of a marching scheme would be inherently unstable, due to the unbounded growth of exp( √ |a| ε x) (in ε). This contrasts with the oscillatory regime, where a basis of the solution would be given by the bounded functions exp(±i √ a ε x). Consequently, we are faced with coupling two different approaches for solving (1.1) in the semi-classical limit: an IVP for a(x) > 0 and a BVP for a(x) < 0.
The goal of this paper is to analyse the numerical coupling of oscillatory and evanescent regimes, using WKB-ansatz functions for both situations. In the oscillatory regime we shall use the marching scheme from [2], and in the evanescent regime a FEM like in Sect. 3 of [19]. In the first case, the key idea is to eliminate analytically the dominant oscillations of the solution to (1.1). The transformed problem then has a much smoother solution, in the sense that the amplitude of the residual oscillations is much smaller than in the original problem -often by many orders of magnitude (in fact by the order ε 2 , cf. Propositions 2.1, 2.2 in [2]). Hence, the new problem can be solved numerically on a coarse grid, still yielding a very accurate approximation. In the evanescent regime, the key idea of the WKB-FEM method is to use WKBansatz functions (of exponential type), rather than the standard polynomials. Finally the strategy is the same as in the oscillatory region: to filter out the boundary layer via well-chosen basis functions. Since WKB-basis functions are asymptotic solutions to (1.1) (as ε → 0), this method is again very accurate on a coarse grid. In this paper we shall provide first the numerical analysis of the WKB-FEM method for the evanescent regime (from scratch), and then the error analysis of the hybrid DDM -building on the convergence results in [2] for the WKB-marching method.
Problems similar to (1.1) or in general that require the numerical integration of highly oscillatory equations play an essential role in a wide range of physical phenomena: e.g. electromagnetic and acoustic scattering (Maxwell and Helmholtz equations in the high frequency regime), wave evolution problems in quantum and plasma physics (Schrödinger equation in the semi-classical regime), and stiff mechanical systems. The application we are interested in here, stems from the electron transport in nanoscale semiconductor devices, like quantum wave-guides [1], resonant tunnelling diodes (RTDs) [6], MOSFETs [21], etc. In a 1D model setting, which is appropriate for RTDs or for the longitudinal dynamics in each transversal mode in MOSFETs, the governing equation is the stationary Schrödinger equation. In an idealized model we assume that the given electrostatic potential V (x) is constant in the left lead x ∈ (−∞, 0], with value V (0), and equally in the right lead x ∈ [1, ∞), with value V (1). Hence the Schrödinger equation can be complemented with open boundary conditions at both ends: This equation describes the state of an electron that is injected with the prescribed energy E > 0 from the right boundary (or lead) into an electronic device (diode, e.g.), modelled on the interval [0, 1]. The corresponding (complex valued) wave function is denoted by ψ(x), where |ψ(x)| 2 is related to the spatial probability density of the electron. Due to the continuous injection of a plane wave function at x = 1, we cannot expect |ψ| 2 to be normalised here (in L 1 (0, 1)). When considering the equivalent problem on R, ψ rather describes a scattering state with ψ ∈ L ∞ (R). The small parameter 0 < ε < 1 is the re-scaled Planck constant. To make the link with (1.1), the coefficient function a(x) is given by a(x) := E − V (x). To allow for an injection at x = 1, we have to require that a(1) > 0, cf. Figure 2. In fact, E > V (x) characterises the oscillatory, classically allowed regime, whereas E < V (x) characterises the evanescent, classically forbidden regime. Fig. 2 sketches a tunnelling structure including both regimes, which are rather different. The boundary conditions in (1.2) are the so called open or transparent boundary conditions, permitting an electron wave to enter or leave the device without reflections [16]. Due to the injected plane wave of electrons, the boundary condition at x = 1 is inhomogeneous.
But at x = 0 it is homogeneous, due to the free outflow of the electron wave.
In the present work we shall not discuss (in detail) situations with turning points, i.e. zeros of a, but rather concentrate on devices with abrupt jumps at the interfaces between oscillatory and evanescent regions. This first step is simpler to treat and will be extended in a subsequent work. In Sect. 5 we shall comment on situations incorporating turning points.
For the solvability of this model, the following simple result holds: Proof For the case of an oscillatory outflow, i.e. a(0) > 0, the proof was provided in Proposition 2.3 of [4]. For an evanescent outflow, i.e. a(0) < 0, the proof is analogous (multiplying the Schrödinger equation byψ, integrating by parts, and taking the imaginary part).

WKB-technique.
Both parts of the hybrid numerical method studied in Sect. 3 will be based on WKB functions. Hence, let us first review the well-known WKB-approximation (cf. [14]) for the singularly perturbed ODE (1.1). In the standard approach, for the oscillatory regime (i.e. a > 0), the WKB-ansatz is inserted in (1.1) and after comparison of the ε p -terms, leads to Truncating the ansatz (1.3) after p = 0, 1, or 2, yields the asymptotic approximation ψ(x) ≈ Cϕ os p (x), with the following oscillatory WKB-functions (of the three lowest orders in ε): Proceeding analogously for the evanescent regime (i.e. a < 0) yields the following evanescent WKB-functions (of the three lowest orders in ε): In the hybrid numerical method analysed in §3 we shall use the first order WKB ansatz functions ϕ ev 1 for the FEM in the evanescent region. And in the oscillatory region we shall use the second order WKB functions ϕ os 2 to transform (1.1) into a smoother problem that can be solved accurately and efficiently on a coarse grid. Since we shall use two different numerical approaches in the two regimes, also the corresponding error orders will be rather different (both with respect to ε and to the grid size h). Hence there is, a-priori, no obvious natural choice for the orders of the two WKB-methods. We choose a first order WKB-method in the evanescent region to keep the complexity of the numerical scheme and the technicalities of its convergence analysis low. Furthermore we choose a second order WKB-method in the oscillatory region such that we can use the results from [2] (without having to redo a first order WKB-analysis). Anyhow, our hybrid method is second order with respect to h. This paper is organised as follows. In Sect. 2 we present and analyse the (nonoverlapping) domain decomposition method for the singularly perturbed ODE (1.1) on the continuous level. Propositions 2.2 and 2.6 establish that this DDM yields the analytical solution in one sweep for cases consisting of two and, resp., three distinct regions. In Sect. 3 we first review the two different numerical WKB-methods for the two distinct regions and establish convergence of the WKB-FEM. Then we prove convergence of the overall hybrid WKB-method (WKB-FEM in the evanescent regime coupled to a WKB-marching scheme for the oscillatory region), with Theorem 3.8 as the main result. In Sect. 4 we illustrate our convergence results on some numerical examples treated with our scheme, including an example with a tunnelling structure. And finally, in Sect. 5 we briefly discuss extensions of our WKB method to coefficient functions with turning points.

Domain decomposition of the Schrödinger boundary value problem
In this section we shall consider the Schrödinger BVP (1.2) with given coefficient functions a(x) = E − V (x) corresponding to two different scenarios -first two coupled regions, then three regions. We confine ourselves here to these examples only for practical reasons: This setup already shows all the interesting features of the BVP, and it can be generalised easily to more regions.

Two coupled regions
We start with the situation illustrated in Fig. 3: It consists of two regimes, an evanescent region with a := E − V < 0 (adjacent to the left boundary) and an oscillatory region with a > 0 (adjacent to the right boundary). Since we exclude turning points here, the function a is assumed to have a jump discontinuity (and a sign change) at the interface x = x d . Moreover, for this section we shall assume: Following the basic idea from [2] we shall solve the BVP (1.2) as two consecutive sub-problems: We start with the evanescent region (0, x d ), where a BVP is solved (for stability reasons, as mentioned in Sect. 1). This is followed by an IVP on the oscillatory region (x d , 1). So we shall proceed in the opposite direction of the injection direction (see Figure 3). Due to Proposition 1.1, the solution to (1.2) is in C 1 [0, 1]. Hence the solutions on these two sub-regions are matched by continuity of ψ and ψ at x = x d .
For the BVP on (0, x d ), the original problem (1.2) provides only one homogeneous Robin boundary condition (BC) at x = 0. Hence, we supply the BVP with an auxiliary, artificial BC at x = x d . Here, both an inhomogeneous Dirichlet or Neumann BC would work from an algorithmic point of view. In order to simplify the numerical anal-ysis in Sect. 3.1 below, we shall use at this point εχ (x d ) = 1 for the auxiliary wave function χ . While this auxiliary value has the correct ε-order (cf. Proposition 2.3 and Lemma 3.4 below), it will in general not be the correct derivative of the global solution ψ. Its correct value will finally be obtained by scaling the auxiliary functions using the remaining inhomogeneous Robin BC at x = 1 (cf. (1.2)). This leads to the following domain decomposition and problem coupling for the auxiliary wave functions χ, ϕ: Step 1 -BVP for χ in region (1): (auxiliary Neumann BC) Step 2 -IVP for ϕ in region (2): Step 3 -Scaling of the auxiliary wave functions: with the scaling parameter α ∈ C defined via , (due to the right BC in (1.2)) .
(2.4) We note that the denominator in this expression for α is non-zero: On the one hand χ and ϕ are both real valued, and on the other hand ϕ(1) and ϕ (1) cannot vanish simultaneously (as otherwise ϕ ≡ 0 would contradict the Neumann BC in (2.1)).
In this whole section we shall only require that a ∈ L ∞ (0, 1). As in Proposition 1.1, a(0) and a(1) are hence not meant as the point values of the function a, but rather as the constant potential in the left and right leads.
The above lemma, whose proof is very easy, shows that the domain decomposition algorithm (2.1)-(2.4) yields a unique function ψ that is piecewise in W 2,∞ and piecewise (on the two regions) a solution to the Schrödinger equation (1.2). In fact, this DDM yields the unique solution of (1.2) as stated in the following proposition: Proposition 2.2 Let Hypothesis a2 be satisfied. Then the function ψ obtained from (2.1)-(2.4) belongs to W 2,∞ (0, 1) and is the unique solution of (1.2) (as guaranteed by Proposition 1.1).
The following result provides the uniform-in-ε boundedness of this solution ψ. It generalizes Theorem 2.2 from [20], which holds only for one purely oscillatory region: Proposition 2.3 Let Hypothesis a2 hold. Moreover, let the potential in the oscillatory region satisfy a ∈ W 1,∞ (x d , 1) and 0 < τ os ≤ a(x) ∀x ∈ (x d , 1). Then, the solution of (1.2) satisfies The simple, but lengthy proof is deferred to the Appendix.

Three coupled regions
In this subsection we consider the Schrödinger BVP (1.2) with a given coefficient function a(x) as illustrated in Fig. 4: It consists of three regimes, two oscillatory regions at the interval boundaries and an evanescent region in the middle. Since we exclude turning points here, a is assumed to have jump discontinuities (and sign changes) at the interfaces x = x c and x = x d . The solution ψ to the BVP (1.2) for such an example is illustrated in Fig. 5 below. Moreover, for this section we shall assume on a(x): Tunnelling structure: while electrons are injected from the right boundary with energy E, the decomposed problem has to be solved from left to right (as an IVP-BVP-IVP). The coefficient function in Following the basic idea from [2] we shall solve the BVP (1.2) as an IVP-BVP-IVP problem in the opposite direction of the injection direction, i.e. starting at x = 0 (see Figure 4). In (1.2), the Robin boundary condition (BC) at x = 0 only fixes the ratio ψ (0) ψ(0) , hence an auxiliary Dirichlet (or Neumann) boundary value has to be invoked here. Its correct value will then be obtained by scaling the final equation using the inhomogeneous Robin BC at x = 1. In contrast to [2], (1.2) includes the evanescent region (2), cf. Fig. 4, which still has to be formulated as a BVP (for stability reasons). This leads to the following domain decomposition and problem coupling for the auxiliary wave functions ζ, χ, ϕ: Step 1 -IVP for ζ in region (1): Step 2 -BVP for χ in region (2): (implies continuity of ψ at x d ) Step 4 -scaling of the auxiliary wave functions: with the scaling parameters α, β ∈ C still to be defined. This procedure can be explained as follows: First we note that the BCs of (2.6) imply εζ (0) + i √ a(0)ζ (0) = 0, just as in the left BC of the BVP (1.2). Hence, the IVP (2.6) coincides on region (1) with the BVP (1.2), except for the auxiliary Dirichlet BC ζ(0) = 1. The true solution of (1.2) satisfies instead ψ(0) = β with some a-priori unknown β ∈ C. Hence, the auxiliary wave function ζ is related to ψ by the scaling ψ [0,x c ] = βζ , as postulated in the first line of (2.9). Clearly, this implies ψ In the above Step 2, the Robin BC allows to carry over this relation to region (2): ψ ψ = χ χ , and the auxiliary wave function χ is related to ψ by the scaling ψ [x c ,x d ] = αχ, with some α ∈ C to be determined. The initial conditions for the auxiliary wave function ϕ in Step 3 imply C 1 -continuity of ψ when using again the scaling So far, the wave function ψ defined in (2.9) neither satisfies continuity at x c nor the right BC from (1.2). Therefore we define the scaling parameters α, β ∈ C via (implies continuity of ψ at x c ) (2.11)

Remark 2.4
The key aspect of the above algorithm is to prescribe in the BVP (2.7) the continuity of ζ ζ to χ χ at x c . Note that this continuity is invariant under the scaling (2.9). Hence it is inherited by ψ ψ , implying (with the continuity of ψ) the required C 1 -continuity of ψ. The simpler alternative to prescribe in (2.7) continuity of ζ to χ would typically be paired with a discontinuity of ζ to χ at x = x c (as a result of solving the BVP). Then, the scaling of (2.9) would lead to an unwanted discontinuity of ψ at x = x c .
Proof Part (i) is straightforward. For (ii), let us first consider the IVP (2.6). Its unique solution ζ has the property: The values ζ(x c ) and ζ (x c ) are linearly independent over the field R. Otherwise, the backward IVP (starting at x c ) would yield "final values" ζ(0) and ζ (0) that are linearly dependent over R, which is in contradiction with the initial condition in (2.6).
To solve the BVP (2.7), let χ 1 , χ 2 be a (real valued) basis of solutions for that Schrödinger equation on (x c , x d ), with Setting χ = c 1 χ 1 + c 2 χ 2 with some c 1 , c 2 ∈ C, the BCs of (2.7) give rise to the following linear equation: and ζ (x c ) are linearly independent over R. Hence, (2.12) is uniquely solvable for c 1 , c 2 , thus providing the unique solution to (2.7). For part (iii) we shall first argue that (2.10) yields a well-defined α ∈ C \ {0}. To this end we shall prove that εϕ (1) Then, (2.9) implies on the one hand But, on the other hand, (2.6) yields Since the current in a stationary quantum model is constant in x, this implies j ≡ 0.
The above lemma shows that the domain decomposition algorithm (2.6)-(2.11) yields a unique function ψ that is piecewise in W 2,∞ and piecewise (on the three regions) a solution to the Schrödinger equation (1.2). Moreover, one has the proposition: Proposition 2.6 Let Hypothesis a3 be satisfied. Then the function ψ obtained from (2.6)-(2.11) belongs to W 2,∞ (0, 1) and is the unique solution of (1.2) (as guaranteed by Proposition 1.1).
Proof The matching conditions in (2.7) and (2.11) imply C 1 -continuity of ψ at x c , and the initial conditions in (2.8) imply C 1 -continuity of ψ at x d . Hence, ψ ∈ C 1 [0, 1], and this proves the claim.

Numerical analysis of the hybrid WKB-method
To keep the presentation simple we shall consider here only the two-zone model of Sect. 2.1. It has a coefficient function a(x) that corresponds to Figure 3. In the following subsections we shall thus study step by step the numerical errors obtained when solving the BVP (2.1) in the evanescent region with a multiscale WKB-FEM and the IVP (2.2) in the oscillatory region with the marching method introduced in [2]. We shall always assume that the phase function x √ |a(y)| dy in the WKB-functions (1.5), (1.6) can be computed exactly, e.g., this holds for piecewise linear a(x). Otherwise, an additional quadrature error of the phase (of the form h γ /ε with γ > 0 denoting the order of the chosen numerical integration) would need to be included in our subsequent analysis. This will be the topic of the follow-up work [3]. That paper will also illustrate how a spectral method for the phase integral allows to drastically reduce this quadrature error in practical computations.

Variational formulation for the evanescent region BVP (2.1)
Let us introduce in this section the variational formulation of the evanescent region problem (2.1) and study the well-posedness of the problem. As pointed out previously, we consider in the current paper situations with an abrupt potential jump, avoiding turning points, such that we shall suppose: Furthermore let in the following 0 < ε < 1 be arbitrary.
We are now searching for a weak solution of (2.1) in the Hilbert space This ε-dependent norm gives rise to the following weighted Sobolev embedding, where the Gagliardo-Nirenberg inequality for bounded domains is used in the first estimate: The variational formulation reads: with the sesquilinear form b : V × V → C and the antilinear form L : The BVP (2.1) is a standard elliptic problem, meaning that the forms b(·, ·) and L(·) are continuous and b(·, ·) is coercive, i.e. there exists a constant C > 0 independent of ε, such that for all χ, θ ∈ V one has The Lax-Milgram theorem implies then for each ε > 0 the existence and uniqueness of a weak solution χ ∈ V of (3.2). We have moreover the following lemma:

Lemma 3.1 Let Hypothesis A be satisfied. Then the weak solution
and satisfies the following estimates, with a constant C > 0 independent on ε

4)
as well as Proof The Lax-Milgram theorem yields immediately which implies (3.4) and, with (3.1), the first inequality of (3.5).
To show the second estimate of (3.5), we observe that where we used the other, just proved estimates.

Convergence analysis for the WKB-FEM method in the evanescent region
The multi-scale WKB-FEM method we shall use for an efficient resolution of the evanescent region problem (2.1) is based on a specific choice of WKB-basis functions from (1.6). In more detail, the Hilbert space V will be approximated by an appropriate finite-dimensional Hilbert space V h ⊂ V, spanned by well chosen basis functions, and the continuous problem (3.2) will be approximated by the following discrete problem: with the WKB-hat functions defined as Here we used the notation Note that both components v n and w n of these (non-standard) hat functions are linear combinations of the evanescent WKBfunctions of first order, i.e. ϕ ev 1 given in (1.6). Hence these hat functions are solutions of our Schrödinger equation up to an error of order O(ε 2 ), i.e.
This peculiarity signifies that the hat functions incorporate already some essential information about the solutions we are searching for, leading to a scheme which will be asymptotically correct in the limit ε → 0, as will be seen later on. For later purposes let us introduce here the differential operator with the function For each ξ ∈ V h one has A ε (ξ ) = 0 in every interval I n := (x n , x n+1 ). But we note that A ε cannot be applied globally on (0, x d ), as functions in V h typically have discontinuous derivatives at the grid points x n . We also remark that, due to Lax-Milgram's theorem, the discrete problem (3.6) admits ∀ε > 0 also a unique solution The aim is now to investigate the error between the exact solution of (3.2), denoted by χ ex , and the solution of the discrete problem (3.6), denoted by χ h . We denote by (3.9) Then, the numerical error can be split as follows where e 1 h corresponds to the interpolation error (consistency) and e 2 h is the stability error. These two error parts shall be now estimated separately.

Consistency error estimate
The goal of this section is to estimate the interpolation error To this end, note that the equation satisfied by e 1 h in I n is The variation of constants method, i.e. writing e 1 h (x) = c 1 (x) w n (x)+c 2 (x) v n (x) in I n leads after some lengthy but straightforward computations (see [20] for the oscillatory case) to the following explicit expressions for the error function with Differentiating the interpolation error function yields with In order to estimate the interpolation error in the V-norm, we shall investigate each of these terms separately. In this study, the behaviour of the following functions is very important: Next we shall use r (y) (V (y)−E) 1/4 ≤ C and the fact that sinh σ n (x) and sinh(−σ n+1 (x)) are both non-negative on I n and, respectively, increasing and decreasing. Then one can show for x ∈ I n : x n+1 x |χ ex (y)| dy, x n+1 x |χ ex (y)| dy + Cε ss (x) x n+1 x |χ ex (y)| dy .
In the above estimates, the constant C depends only on our data a(x) and E, but not on ε and χ ex . Using σ n (x) − σ n+1 (x) = γ n we easily find With the asymptotic behaviour With Lemma 3.1 this permits to prove the following lemma:

Stability error estimate
Let us now come to the V-estimates of the stability error For this study, we remark that Using now the coercivity and the boundedness of the sesquilinear form b, yields with a constant C > 0 independent of ε, that but this estimate can be further improved for 0 < ε 1. For this, let us study the right hand side of (3.12) in more detail. As e 2 h ∈ V h we have that A ε (e 2 h ) = 0 on every interval I n , implying thus were we used e 1 h (x n ) = 0 for the integration by parts. Thus one obtains from (3.3), (3.12): and in particular e 2 h L 2 ≤ Cε 2 e 1 h L 2 . Using also the Sobolev embedding (3.1), this implies Finally we shall now derive an L ∞ -bound on (e 2 h ) , using the bound on To bound ζ n L ∞ , one refers to (3.8) and observes that 0 ≤ α n (x) ≤ 1 as well as And analogous estimates hold for β n (x). From (3.13) we then obtain (3.14)

Convergence results for the WKB-FEM method
To summarize, we shall put together both error contributions (from Lemma 3.2 and from the stability estimates (3.13), (3.14)) in the following theorem. It turns out that the consistency error is dominant here:

Theorem 3.3 (Convergence WKB-FEM)
Let Hypothesis A be satisfied. Then the following estimates hold for the numerical error between the exact solution χ ex ∈ V of (3.2) and the numerical solution χ h ∈ V h of (3.6):

Vectorial IVP for the oscillatory region
In this subsection we shall first investigate the IVP (2.2) on the continuous level, on the interval (x d , 1). Following [2], we shall rewrite (2.2) in vectorial form, done via a non-standard transformation that is appropriate for the numerical WKB-marching method. For the subsequent analysis let us make the following assumptions on the potential: 1] and E > 0 satisfy the bounds Moreover let 0 < ε ≤ ε 0 be arbitrary, with some ε 0 such that 0 < ε 0 < ε 1 := min 1, min In this definition, β + denotes the non-negative part of β. Hence β + (x) −1/2 may take the value ∞. We note that the above restriction on ε guarantees that the phase function of ϕ os 2 (x) (cf. (1.5)) is strictly increasing. Moreover, the resulting positivity of the function √ a − ε 2 β will be crucial for the WKB-marching method in Sect. 3.4. Following [2] it is convenient to pass from the second-order differential equation to a system of first-order, introducing the following vector notation for the wave function ϕ(x) on [x d , 1]: . (3.15) The norm of U is equivalent to the norm of the vector (ϕ, εϕ ) . Indeed, the transformation matrix between these two vectors reads i.e. U (x) = A(x) ϕ εϕ , (3.16) where the matrix A and its inverse are bounded, uniformly w.r.t. x and ε, due to Hypothesis B. Let ϕ ex ∈ W 2,∞ (x d , 1) be the exact solution of (2.2) as guaranteed by Lemma 2.1. In the above vector notation it will be denoted by U ex (x) or simply U (x), and is solution to the system with the two matrices Here, β = − 1 2a 1/4 (a −1/4 ) which was already defined in (1.4), and the matrix element a(x d +) of A(x d +) denotes the right-sided limit of a at the jump discontinuity x d . We also use the analogous notation for a (x d +).
In the sequel we shall need an a-priori estimate on this solution. The upper bound was already given in §2.1 of [2]. But for the scaling Step 3 we shall also need an εuniform lower bound on the solution. Here and in the sequel we shall use the notation U 2 := |u 1 | 2 + |u 2 | 2 .

Lemma 3.4 Let Hypothesis B hold. Then, the ODE (3.17) admits a unique solution
Thus, there exist constants C 3 , C 4 > 0 independent on ε such that Proof For the norm U 2 we compute for (3.17): This implies (3.18). The estimate (3.19) is now a simple consequence of (3.18), presupposing that one proves some ε-independent bounds on the initial condition U ex (x d ) , or equivalently (χ ex (x d ), 1) . The latter norm is clearly bounded below by 1, and it is also bounded above due to the a-priori estimate on χ ex C[0,x d ] from Lemma 3.1.

Review of the WKB-marching method for the oscillatory region
In this subsection we shall first review the WKB-marching method for solving the IVP (2.2) (or, equivalently, (3.17)). Then we recall its error estimates from [2]. Following [2] this method consists of two parts, first an analytic transformation of (2.2) or (3.17) into a less oscillatory problem, and second the discretization of the smooth problem on a coarse grid in an ε-uniform manner. As shown in [2], the analytic WKB-transformation reviewed here is related to using oscillatory WKB-functions of second order, ϕ os 2 (x). Part 1 -analytic transformation: The starting point is the vectorial IVP (3.17).
The vector function U ∈ C 2 is then transformed to the new unknown Z ∈ C 2 by (3.21) with the matrices , and the (real valued) phase function Note that this is precisely the phase in the second order WKB-approximation ϕ os 2 (x) from (1.5). In (3.21), the unitary matrix P is chosen to diagonalise A 0 , i.e. the dominant part of (3.17) as ε → 0. Next, the diagonal matrix exp(− i ε ε ) eliminates the leading oscillations.
This change of unknown leads to the smooth ODE-system (3.23) Here, the 2 × 2-matrix function is bounded independently of ε. It is off-diagonal, with the entries This finishes the analytical transformation, and the goal of the second part is to provide an ε-uniform discretisation of (3.23) that is second order w.r.t. the mesh size. With the initial condition Z N := P U N ∈ C 2 and U N := U ex (x d ) ∈ C 2 given, the marching scheme reads as follows (see [2]): Before listing the two matrices A 1 n , A 2 n we give a short motivation for their derivation: This rather non-standard method was constructed such that it is at the same time second order convergent (w.r.t. h) and uniform w.r.t. ε. For each marching step we first use a second order Picard approximation of (3.23), which leads to ε-oscillatory integrals in x (due to the oscillations of N ε (x)). These integrals involve the two small parameters ε and h, and are approximated using variants of the asymptotic method [12].
In (3.24) the 2 × 2-matrices are given by Here we used the notation and the discrete phase increments Remark that for notational reasons we omitted in the aforementioned description of the scheme the ε-index. Furthermore we assumed that the two functions φ and β (the latter involving the derivatives a , a ) are explicitly "available". Alternatively, φ, a and a could be approximated numerically. But, for simplicity, we shall not include such errors in the subsequent error analysis. Proof Let us start by analysing the propagation matrix B n := I + A 1 n + A 2 n ∈ C 2×2 of the vector Z n as defined in (3.24). A straightforward computation reveals its symmetry (which was also used in the proof of Proposition 3.3, [2]): This symmetry carries over to the matrix With this notation, the propagation matrix for the vector U n reads (cf. (3.24), (3.25)): where we used This shows that U n ∈ R 2 . Coming now to the bounds of U n , a simple Taylor expansion for the matrices in (3.24) yields A 1 n ≤ Cεh n , A 2 n ≤ Cε 3 h n , and hence with some constant C 7 > 0: For the lower bound we chooseε 0 := min{ε 0 , With the elementary estimate 1 1−y ≤ 4 y for y ∈ [0, 0.5] we then obtain This allows to estimate the backwards propagation Z n = ( and the lower bound on U n follows as before.
Due to the above lemma we have to restrict the range of admissible ε-values: Hypothesis B' Let the assumptions of Hypothesis B hold, but with ε 0 replaced byε 0 from Lemma 3.5.

Error and stability estimates for the WKB-marching method
In this subsection we recall the main Theorem 3.1 of [2], providing error bounds for the marching method (3.24)-(3.25), used for solving the IVP (2.2) or, equivalently, (3.17).

27)
with C independent of n, h, and ε. Here, γ > 0 is the order of the chosen numerical integration method for computing the phase integral for the back-transformation (3.25).
We remark that the term h γ ε of (3.27) may or may not be present in real computations, depending on the chosen coefficient function a(x). If a(x) is piecewise linear or piecewise quadratic, e.g., the phase integral φ ε (x) can be computed analytically. Hence, this term would not appear in such cases. In the numerical tests performed in Sect. 4 below, we shall only consider such examples of exactly computable phase functions and shall hence not include this error term in the error analysis of Sect. 3.5.

Convergence results for the overall hybrid WKB method
In this section we shall combine the error analysis of the previous two sections and adapt it to the algorithm for coupling two regions. To this end we have to fix the numerical analogues of the continuous coupling conditions in (2.2). First we shall (of course) use ϕ h (x d ) := χ h (x d ). But for the initial condition of the derivative there are two options, namely εϕ h (x d ) := εχ h (x d ) or εϕ h (x d ) := 1 (taken from the exact value in (2.1)). We shall use the second option for the following reasons: On the one hand it avoids the numerical error of χ h (x d ), where we recall that χ h is discontinuous (and hence less accurate) at the grid points. And on the other hand, this choice will facilitate the a-priori estimate needed for the scaling in Step 3.
Since the numerically used initial data χ h (x d ) deviates from its exact value χ ex (x d ), this gives rise to an additional error component to be considered: Let thusÛ ex (x) denote the exact solution to the ODE (3.17), but with the following perturbed initial condition:Û Using the a-priori estimate (3.18) leads to the error where we used the notation e h = χ ex − χ h . The following convergence analysis of the hybrid method uses several different solution functions (exact, numerical, etc.). To keep the notation straight we summarize it in the following table, both for the evanescent and oscillatory regions. The superscript ( ) signifies that we refer to both the function and its first derivative. Numerical solution vector after scaling in Step 3, with numerical IC For clarity, we summarize here also the numerical analogue of the three steps in (2.1)-(2.4), referring to the two regions in Fig. 3: Step 1 -WKB-FEM for χ h in region (1): Step 2 -WKB-marching method for ϕ h in region (2): As initial condition for the marching scheme we use U N :=Û ex (x d ) ∈ R 2 given by (3.26). Applying the scheme (3.24)-(3.25) iteratively we compute the vectors U n = (u 1 n , u 2 n ) ∈ R 2 ; n = N + 1, . . . , M.
Step 3 -Scaling of the auxiliary wave functions χ h , ϕ h : with the scaling parameterα ∈ C defined in analogy to (2.4): The statement (3.32) reveals that our final numerical solution ψ h is continuous in the evanescent region, but discrete in the oscillatory region. Note also that the relation between U = (u 1 , u 2 ) and (ϕ, εϕ ) , given by (3.16), provides a connection between the two scaling functions α andα, i.e.
Let us finally also recall that the solution χ h as well as the vector U n = 0 of Step 2 are real valued. The final (numerical) solution ψ h only becomes complex valued due to the multiplication byα in Step 3. Note that the denominator of (3.33) cannot vanish for U ∈ R 2 \{0}, which makes the scaling well defined. This mapα satisfies moreover the following simple properties: Lemma 3.7 For each fixed δ > 0, the (ε-dependent) mapα : U ∈ R 2 \ B δ (0) → C is Lipschitz continuous with some constant L > 0 and bounded by some constant C 8 . Both constants can be chosen uniformly w.r.t. 0 < ε ≤ ε 0 and are dependent on δ.
The following error analysis of the hybrid scheme is the main result of this paper.  with the notationẽ h,n := ψ ex (x n ) − ψ h,n andẽ h,n := ψ ex (x n ) − ψ h,n . c) For the overall hybrid method one has then, over [0,1], the estimates We continue with estimating the error propagation due to the non-linear scaling in Step 3. Due to Lemma 3.7, the mapα : U ∈ R 2 \ B min{C 3 ,C 5 } (0) → C is Lipschitz continuous with some constant L > 0 and bounded by some constant C 8 . Both of these constants can be chosen independent of 0 < ε ≤ ε 0 , as the choice δ := min{C 3 , C 5 } for the domain ofα uses the lower bounds on U ex from (3.19) and on U n from Lemma 3.5.
Here it is crucial that both the exact solution vector U ex (1) This yields where we used (3.37) in the last step. Using now Lemma 3.1 to estimate χ ex ( ) (x) and Theorem 3.3 for e h ( ) (x) yields the four error estimates of (3.34).

Part b)
For the oscillatory region [x d , 1] we estimate the difference between the exact solution and the numerical solution (both after scaling) at the grid points x n ; n = N , . . . , M. Here it is again more convenient to use the vector notation from (3.15), where we introduce the notationsŨ ex (x) :=αU ex (x) andŨ n :=αU n for the exact and, respectively, numerical solution after scaling in Step 3: where we used (3.19) in the penultimate line, and (3.37) twice in the last line. Using the norm equivalence then yields the error estimate (3.35).
Part c) is just a combination of the previous two parts.

Numerical tests of the hybrid WKB method
The aim of this section is to present numerical results obtained with the WKB-coupling scheme introduced in Sect. 2 and to compare these results with the error analysis established in Sect. 3.5. In particular, we present the results for 3 zones (oscillatingevanescent-oscillating, cf. Sect. 2.2) corresponding to the passage or flow of electrons through a tunnelling structure (see Fig. 5), with a piecewise linear and, respectively, piecewise quadratic potential V (x), chosen such that the phase φ ε (x) is explicitly calculable. The reason for such a choice is to avoid having to care about the h γ ε -error term in (3.27), yielding hence an asymptotically correct scheme for fixed h > 0 and ε → 0.
In Fig. 6 we plotted the numerical errors of the coupling method associated to the wave function ψ (left figure) and to its derivative εψ (right figure), as functions of the mesh size h (in log − log scale) and for three different ε-values. In the oscillating regions we chose the second order method (3.24)-(3.25) and in the evanescent region the FEM (3.6). The plotted errors are the L ∞ -errors between the numerical solution on the whole interval [0, 1] and a reference solution, computed with the same scheme but on a finer grid of 2 18 points. It can be observed that the slopes in these two plots are approximatively one (for h 3 · 10 −5 ) and improving to 1.5 for smaller values of Summarizing, the error of ψ shows to be of order O(min{h 2 , √ εh 3/2 , εh, ε 3/2 √ h}), which corresponds exactly to the estimates given in Theorem 3.8(c). The error of εψ shows to be of order O(min{ √ εh 3/2 , εh, ε 3/2 √ h}), which is even slightly better than the estimate from Theorem 3.8(c) (in the sense of including also an O( √ εh 3/2 )behaviour).
We mention that the obtained numerical errors are mainly those introduced by the WKB-FEM of the evanescent region. Indeed, in this evanescent region, the numerical error of the WKB-FEM is larger than the one obtained from the second order WKB marching method of the oscillating region (compare the estimates in the Theorems 3.3 and 3.6).
Example 2 Next we consider a piecewise quadratic potential given by with x c = 0.5, x d = 0.5 + 2 −5 = 0.53125 and Before turning to the error plots we consider the condition number associated to solving the discrete variational problem (3.6) in the (intermediate) evanescent region. In Fig. 7 we plot this condition number as a function of h, for three different values of We remark that this behaviour is not a problem in practice: For large ε, the solution ψ ex is not highly oscillatory and hence does not need a high spatial resolution. For small values of ε, even a fine resolution would only lead to moderate condition numbers. Indeed, one observes a decrease of the condition number when ε gets smaller. This important feature is somehow related to the asymptotic-preserving property of the scheme. Large condition numbers signify that the errors of numerical experiments also include significant contributions stemming from round-of errors and their accumulation. While the method-error (as estimated in Theorem 3.8(c)) decreases with decreasing h, the round-of errors could then increase in some situations, due to the increasing condition number. These arguments may lead to the idea that one cannot trust too much the reference solution in Fig. 6, computed with 2 18 points. In order to verify this suspicion, we decided to plot in the case of a piecewise quadratic potential two types of error curves.
To be more precise, in Fig. 8 we show the numerical errors of the wave function ψ (left figure) and its derivative εψ (right figure), as functions of the mesh size h (in log − log scale) and for four different ε-values. The dashed lines correspond (as in Example 1) to the L ∞ -error between the numerical solution on the whole interval [0, 1] and a reference solution, computed with the same scheme but on a finer grid (here h = 2 −19 ) whereas the solid lines correspond to the incremental error when iteratively doubling the grid size, i.e. ||ψ h j − ψ h j−1 || ∞ with h j = 2h j−1 . For a first order method, the former error is about twice as large as the latter (incremental) error. This can be understood from the geometric series of the incremental errors, since the summands then differ by a factor of about 2. In Fig. 8 this difference is clearly visible for the solid red curves, pertaining to ε = 10 −2 , and the corresponding dashed error curves (for large h). The minimum of the incremental error (as a function of h) indicates the onset of significant round-of errors when reducing h. In Fig. 8 this is best visible for the solid blue and red curves, pertaining to ε = 10 −1 , 10 −2 .
Furthermore remark that for ε = 10 −1 , 10 −2 , 10 −3 and h 3 · 10 −5 , the error slopes are approximately one-just like in Example 1. For smaller values of h the error then gets polluted by round-of errors. For ε = 10 −4 the shown errors seem to be mostly due to round-of errors. They again increase for h 3 · 10 −5 .

Turning points
A turning point of the Schrödinger equation (1.1) is defined as a zero of the given coefficient function a(x). Accordingly one also speaks about the order of a turning point. We first remark that both error analyses, in Sect. 3.2 for the WKB-FEM and in [2] for the WKB-marching method are not valid for turning points. Therefore, we assumed in Hypothesis A and B that a(x) is bounded away from zero. Furthermore, in the convergence Theorems 3.3 and 3.6 we did not keep track how the leading constant C grows with τ ev , τ os → 0. However, the failure of both WKB methods when approaching turning points also appeared in our numerical experiments (not included here). The paper [17] considers a matrix generalization of our equation (1.1), but only for the oscillatory case. Their coefficient matrix A(x) (generalizing our a(x)) is there assumed to be symmetric positive definite, satisfying the uniform lower bound A(x) ≥ δ 2 > 0. The proof of their Theorem 6.1 shows that their L ∞ -error bounds would grow like O(δ −2 ). This numerical failure can be understood easily: The WKB-ansatz (1.3)-(1.6) is not valid at turning points. In fact, close to a turning point of first order, solutions to (1.1) are neither exponential nor oscillatory, but in a transition layer of thickness O(ε 2/3 ) they behave rather like Airy functions: Indeed, for a(x) := − x, a solution basis for (1.1) is given by Ai(ε −2/3 x), Bi(ε −2/3 x), where Ai and Bi denote the Airy functions of first and second kinds.
At a turning point the solution to (1.1) clearly satisfies ψ = 0. This motivated to use linear FEM-ansatz functions in the numerical cell containing a turning point (cf. §3.2.2 in [19]).
The quest for an appropriate replacement of the WKB-ansatz close to turning points has a long history in asymptotic analysis, starting with Langer [15]: For general coefficients a(x) with a zero, he found an asymptotic approximation for the solution of (1.1) that is valid uniformly in x, including the turning point. For a first order turning point his approximation is a composite function involving Airy functions and the phase function (like σ (x) defined in (3.8)). For details on first order turning points we refer to §4.3 in [9], and to §7.3 in [18] for higher order turning points.
The above mentioned approximation formulas of Langer have also been used for numerical computations, mostly for Schrödinger eigenvalue problems [7,8,22]. In the physics literature, this strategy is frequently called Modified Airy function (MAF) technique. It relies on evaluating the explicit formulas of approximate solutions, but it has not been the starting point of constructing a (convergent) numerical method. In a follow-up paper we shall use Langer's approximation functions as ansatz-functions for an ε-uniform numerical method that should also cover turning points.

Conclusion
This paper is concerned with a 1D Schrödinger scattering problem in the semi-classical limit, with the inflow given by plane waves. The injection energy and potential are given such that the problem involves both oscillatory and evanescent regions. For the continuous boundary value problem we presented a new, non-overlapping domain decomposition method that separates the oscillatory and evanescent subproblems. The former are treated as IVPs, and the latter as BVPs. Key issues of this approach are the appropriate interface conditions and the final scaling of the solution function. We proved that the domain decomposition method yields the exact solution in a single sweep, performed in the opposite direction of the wave injection.
The hybrid numerical discretization is based on WKB-methods in both types of regions: a WKB-FEM for evanescent regions [20], and a WKB-marching method for oscillatory regions [2]. The objective of these WKB-methods is to provide an accurate solution -even on coarse grids and independently of ε. Hence, they are asymptotic preserving. For the first time we present an error analysis for the WKB-FEM method. Together with the analysis of the WKB-marching method from [2], this constitutes the key ingredient for our complete convergence analysis of the hybrid WKB-method. Finally, these error bounds are illustrated and verified in numerical experiments. Since the current in a stationary quantum model is constant in x, this implies |r | = 1. Then (7.1) implies With this bound for the initial condition at x = 1 we now consider the IVP (1.1) on [x d , 1]. Then, Theorem 2.2 from [20] yields the asserted estimate (2.5) on the oscillatory region.
Step 2: For the evanescent region [0, x d ] we consider the scaled, real valued solution χ of the BVP (2.1). With an elementary argument we first show that χ has no zero in [0, x d ]: Assuming the opposite, let χ(x 0 ) = 0, which implies χ (x 0 ) = 0 (as otherwise χ ≡ 0). Then χ is convex on one side of x 0 and concave on the other side, with sgn(χ ) = sgn(χ ) due to a [0,x d ) < 0. But then χ(0) and χ (0) would have opposite signs, contradicting the left BC in (2.1). So we conclude that χ does not change signs in [0, x d ].
After scaling this auxiliary function, we find that also |ψ| and |εψ | are increasing on [0, x d ]. Therefore the uniform bound (2.5) carries over from x = x d to the evanescent region [0, x d ].