Stationary Schr\"odinger equation in the semi-classical limit: WKB-based scheme coupled to a turning point

This paper is concerned with the efficient numerical treatment of 1D stationary Schr\"odinger equations in the semi-classical limit when including a turning point of first order. For the considered scattering problems we show that the wave function asymptotically blows up at the turning point as the scaled Planck constant $\varepsilon\to0$, which is a key challenge for the analysis. Assuming that the given potential is linear or quadratic in a small neighborhood of the turning point, the problem is analytically solvable on that subinterval in terms of Airy or parabolic cylinder functions, respectively. Away from the turning point, the analytical solution is coupled to a numerical solution that is based on a WKB-marching method -- using a coarse grid even for highly oscillatory solutions. We provide an error analysis for the hybrid analytic-numerical problem up to the turning point (where the solution is asymptotically unbounded) and illustrate it in numerical experiments: If the phase of the problem is explicitly computable, the hybrid scheme is asymptotically correct w.r.t.\ $\varepsilon$. If the phase is obtained with a quadrature rule of, e.g., order 4, then the spatial grid size has the limitation $h={\cal O}(\varepsilon^{7/12})$ which is slightly worse than the $h={\cal O}(\varepsilon^{1/2})$ restriction in the case without a turning point.

1. Introduction. This paper is concerned with the numerical treatment of the stationary, one-dimensional Schrödinger equation (1.1) ε 2 ψ (x) + a(x)ψ(x) = 0 , x ∈ R in the semi-classical limit ε → 0. Here, 0 < ε 1 is the rescaled Planck constant (ε := √ 2m ), ψ(x) the (possibly complex valued) Schrödinger wavefunction, and the (given) real valued coefficient function a(x) is related to the potential. For a(x) > 0 and ε "small", the solution is highly oscillatory, and efficient numerical schemes have been developed, e.g., in [ABN11,LJL05]. The novel feature of this work is to include one turning point of first order at x = 0 (i.e. a(0) = 0, a (0) > 0). Then, x < 0 represents the evanescent region where the solution ψ decreases exponentially, possibly including a pronounced boundary layer. x > 0 is the oscillatory region, where the solution exhibits rapid oscillations with (local) wave length λ(x) = 2πε √ a(x) . Hence, the situation at hand is a classical multi-scale problem. Standard numerical methods (e.g., [IB95,IB97]) for (1.1) -particularly in the highly oscillatory regime-are costly and inefficient as they would require to resolve the oscillations by choosing a spatial grid with step size h = O(ε). By contrast, we are aiming here at a numerical method on a coarse grid with h > λ, while still recovering the fine structures of the solution. Our strategy is built upon the following works: For the purely oscillatory case a(x) ≥ τ 1 > 0 in the semi-classical regime, WKBbased marching methods (named after the physicists Wentzel, Kramers, Brillouin) physics). Our interest in this problem is motivated by the electron transport in nanoscale semiconductor devices, which will determine the details of our set-up. 1D models are of course idealizations, but quite appropriate, e.g., for resonant tunneling diodes [BP06]. In this application ψ represents the quantum mechanical wave function. The derived macroscopic quantities of interest to practitioners are the particle density n(x) := |ψ(x)| 2 and the current density j(x) := ε (ψ(x)ψ (x)), see [ABN11] for more details. While ψ is highly oscillatory, n and j are not; but they can only be obtained from the wave function.
We consider the internal domain x ∈ [0, 1], which corresponds to the semiconductor device. Moreover, we assume that electrons are injected from the right boundary (or lead) with the prescribed energy E. The coefficient function a in (1.1) is then given by a(x) := E − V (x) where V (x) is the (prescribed) electrostatic potential of the problem. In reality, electrons are injected into a device as a statistical mixture with continuous energies E ≥ E 0 [FG97]. Therefore, given a potential V , a whole energy interval will give rise to turning points 1 , i.e. zeros of a within the interval [0, 1]. To simplify the presentation we shall assume that the only turning point is located at x = 0, and we shall consider only one such injection energy E. Specifically, we shall make the following assumptions for the given potential V (see Fig. 1).
Assumption 2.1. a) V ∈ C(R; R). Potential jumps are excluded here only for simplicity. Without difficulty, they could be included inside the interior domain (0, 1) by restarting the IVP (from Section 3) at jump points. b) Let V (x) < V (0) for x ∈ (0, 1] (to exclude further turning points besides of x = 0). c) The potential in the left exterior domain (i.e. x < 0), and also in a (small) neighborhood of x = 0, is linear (for simplicity of the hybrid problem). We also assume (w.l.o.g.) that it has slope −1. More precisely we assume that ∃ x 1 ∈ (0, 1) such that a(x) = x for x ≤ x 1 . Hence, (1.1) is the (scaled) Airy equation for x ≤ x 1 . d) The potential in the right exterior domain (i.e. x > 1) is constant with value V (1).
Hence, this scattering problem is oscillatory for x > 0, evanescent for x < 0, and it has a turning point of order 1 at x = 0. Since V (x) > E for x < 0, the injected wave function is fully reflected. Due to a) and b) ∃ τ 1 > 0 such that τ 1 ≤ a(x) on [x 1 , 1], which is an important assumption for the WKB-marching method from [ABN11].
In a realistic device model, it would of course be appropriate to assume that the potential V is constant also in the left lead, i.e. on (−∞, x 2 ] with some x 2 < 0. If V (x 2 ) > E, the injected wave would still be fully reflected, leading to a situation that is qualitatively very similar to the present case. In order to simplify the subsequent proofs, we shall stick here to Assumption 2.1 c). The extension to a constant potential on (−∞, x 2 ] will be discussed in a follow-up work. For each fixed ε, the wave function ψ is C 2 and we require the scattering solution to be bounded w.r.t. x ∈ R. Since the potential grows linearly for x < 0, ψ(x) has to decay to 0 as x → −∞. Hence, ψ is a scaled Airy function Ai for x ≤ x 1 : E 0 explicit sol. x 1 numerical solution 1 expl. sol. with some c 0 ∈ C to be determined. Moreover ψ is a superposition of two plane waves (incoming and outgoing) for x ≥ 1: with a 1 := a(1). This whole-space problem (with x ∈ R) can be written as an equivalent BVP for ψ on the interval [x 1 , 1] by using transparent boundary conditions (BCs) that correspond to the C 1 -continuity of the matched whole-space solution. The inhomogeneous transparent BC at x = 1 is well known from [LK90, ABN11] (see (2.4)) and ensures that there is no reflection induced by the BC for an incident wave coming from the right. C 1 -matching of ψ with ψ − at x 1 reads: with some c 0 (ε) ∈ C \ {0}. Eliminating the (so far) unknown constant c 0 , the last two conditions are combined into a Robin BC. In summary this yields the following BVP: Here we already used the assumption that the incident wave has amplitude 1, and more precisely that c 2 = e i √ a 1 ε . Since the wave is fully reflected, we have |c 1 | = 1 for the reflection coefficient c 1 . Hence, the wave ψ + from (2.2) has the maximum amplitude 2 on x ≥ 1. The plots in Figures 2 and 7 illustrate this.
Proof. This proof is analogous to Proposition 2.3 of [BDM97] and Proposition 1.1 of [AN18]: multiplying the Schrödinger equation byψ, integrating by parts, and taking the imaginary part.
Notation and assumptions: Now we recall some notation and assumptions needed to apply the WKB-marching method from [ABN11] to the BVP (2.4). The well-known WKB-approximation (cf. [LL85], §15 of [Hal13]), for the oscillatory regime where a(x) ≥ τ 1 > 0, is based on inserting the asymptotic power series ansatz into the equation (1.1), and comparing ε-powers to successively obtain the functions φ p (x). Truncating the sum in the exponential after p = 2 leads to the 2 nd order asymptotic WKB-approximation for the oscillatory regime with the phase In the WKB-marching method from [ABN11] this 2 nd order WKB-approximation is used to transform the equation (1.1) to a smoother problem that is then numerically solved on a coarse grid, accurately and efficiently. This is done by reformulating the BVP (2.4) into an IVP using the boundary condition at x 1 and then scaling the numerical approximation to this IVP (obtained by the scheme we will recall in Section 5.1) to also satisfy the boundary condition at x = 1. We need to make an assumption to assure the feasibility of the WKB-marching method: Assumption 2.3. Let a ∈ C 5 [x 1 , 1] be real valued and satisfy the following bounds where β + denotes the non-negative part of β.
Mind that this assumption on ε guarantees that the phase φ(x) for the 2 nd order WKBapproximation is strictly increasing since the integrand √ a − ε 2 β is then positive.
3. Analytical problem: reformulation as IVP. For the numerical solution of (2.4) we want to apply the WKB-marching method from [ABN11] on the interval [x 1 , 1]. To this end we need to reformulate the BVP (2.4) as an IVP, whose solution ψ will be scaled afterwards to satisfy the BCs in (2.4). Note that the left BC at x 1 in (2.4) is invariant under scalings.
Step 1: Since (2.4) only includes one condition at x 1 , it is necessary to prescribe an additional, auxiliary initial value forψ(x 1 ). The condition on εψ (x 1 ) then follows from the Robin BC at x 1 in (2.4). On the one hand the two ICs should have the structure of (2.3) (with an appropriate choice of c 0 ), and on the other hand the scaling constant c 0 should be of order 2 Θ(ε − 1 6 ), such that the initial condition vector (ψ(x 1 ), εψ (x 1 )) is ε-uniformly bounded above and below. In fact, c 0 = ε − 1 6 in (2.3) yields where the constants c, C > 0 are independent of ε. This can be verified using the asymptotic expansions for Ai(−z) and Ai (−z) from (A.4): We get the asymptotic representations with the argument z = x ε 2/3 for some (fixed) x > 0 as ε → 0: The ε-uniform lower bound on the IC can be found, as cos and sin never vanish simultaneously. Hence, this scaling gives a natural balance of ψ − and εψ − . We shall thus consider the IVP Here and in the sequel we use the notationψ to refer to the solution of this IVP.
Step 2: Next the solutionψ of this IVP is scaled as in order to satisfy the BC at x = 1. Note that this scaling preserves the BC at x 1 , and thus, ψ is a solution to the BVP (2.4). From here on we use the notation ψ to refer to the solution of BVP (2.4). This scaling is also applied to the extension of the solution to [0, x 1 ] as ψ(x) = α ψ − (x), with the choice c 0 = ε − 1 6 in (2.1). This equivalence of the BVP to an IVP with a-posteriori scaling was already used in [ABN11,§2] and [AN18, Prop. 2.2] for closely related problems. The vector valued system: Following [ABN11] it is convenient to reformulate the second order differential equation (1.1) as a system of first order. This is done in the following non-standard way: Instead of the vector (ψ(x), εψ (x)) we shall use , with the transformation matrix Under Assumption 2.3 (i.e. a(x) is bounded away from zero), the transformation matrix A(x) and its inverse are uniformly bounded w.r.t. x ∈ [x 1 , 1] and ε. Hence, the norms of the two vectorsŴ (x) and (ψ(x), εψ (x)) are equivalent, uniformly in ε. After the transformation (3.6), the IVP (3.3) reads with the two matrices In order to show (in Section 5) that the WKB-marching method from [ABN11] applied to (3.3) yields a uniformly accurate scheme for the BVP (2.4), we shall need ε-uniform boundedness of the scaling factor α from (3.5). This can be inferred from a uniform lower bound on (ψ, εψ ) , which we establish similarly to [AN18, Lemma 3.4]: Lemma 3.1. Let a(x) ∈ C 2 [x 1 , 1] and a(x) ≥ τ 1 > 0. Letψ(x) be the solution to the IVP (3.3), then (ψ(x), εψ (x)) is uniformly bounded above and below, i.e. (3.10) where the constants C 1 , . . . , C 4 > 0 are independent of 0 < ε < ε 0 .
A solutionŴ to the analytical IVP (3.8) needs to be scaled such that, after transforming back via A(x) −1 to (ψ, εψ ) , it fits both boundary conditions in (2.4). In analogy to (3.5) this is done via the scaling parameterα ∈ C defined as (3.12)α(Ŵ (1)) : Now we can write the exact solution ψ(x) to the BVP (2.4) extended to the region [0, 1] as W , as well asψ, are real-valued, and ψ only becomes complex-valued due to the scaling withα(= α). SinceŴ ∈ R 2 , the denominator inα cannot vanish, except for the trivial solutionŴ ≡ 0. Therefore the scaling byα is well-defined and one can show the following properties.
4. Asymptotic blow-up at the turning point. The goal of this work is to construct an ε-uniformly accurate numerical scheme for (2.4). Since it incorporates (implicitly) the turning point at x = 0, it shall be a generalization of [ABN11,AN18].
A key ingredient of the numerical analysis in these papers was the uniform boundedness of the solution ψ w.r.t. ε. But when including a turning point, this does not hold any more, which is a main challenge of the situation at hand. At the turning point, solutions to the BVP (2.4) exhibit blow-up behavior as ε → 0, i.e. |ψ(0)| = Θ(ε − 1 6 ). This is shown in the following explicitly solvable example and the proposition that follows it.
Proof. For readability of this proof, we omit the index ε in ψ ε andψ ε . As discussed in Section 3, the solution ψ to the BVP (2.4) can be obtained by scaling the solution ψ to the IVP (3.3) with the constant α from (3.5), i.e. ψ = αψ. With some (fixed) x 0 ∈ (0, x 1 ) we make a case distinction: Region x 0 ≤ x ≤ 1: For x ≥ x 0 we consider the IVP (3.8) first with a generic initial condition at x 0 . For such x we have a(x) ≥ τ 3 > 0, and therefore the transformation matrices A(x) and its inverse are uniformly bounded w.r.t. x and ε. Denoting by W the vector valued solution to this IVP, this matrix bound implies equivalence (uniform in ε) of the norms of the two vectors W (x) and (ψ(x), εψ (x)) . Note that this equivalence would not hold for the choice x 0 = 0. Analogously as in the proof of Lemma 3.1 we obtain Step 1: First we consider an auxiliary problem that corresponds to the "pure" Airy function solution: LetW be the solution of the IVP (3.8) on x ≥ x 0 with ε = 1, a(x) = x and the initial conditioñ Then the estimates (4.2) with with the constants defined as |β(y)|dy .
Remark 4.3. The proof of the above proposition illustrates the reason for the asymptotic blow-up of |ψ ε (0)|: Essentially, it stems from the x −1/4 -decay of the flipped Airy function Close to the turning point of first order, ψ ε behaves like the scaled Airy function ε −1/6 Ai(−x ε −2/3 ). In the scattering model of Section 2-5 this even holds exactly, and in Section 6 this will hold approximately. This ε-scaling of the x variable compresses this Airy function decay to the (small) interval [0, x 1 ]. At the fixed point x 1 , Ai(−x 1 ε −2/3 ) is proportional to ε 1/6 , which we compensated by the scaling ε −1/6 yielding an ε-uniformly bounded initial conditionŴ (x 1 ). This ε-uniformity is then not affected any more on the subsequent interval [x 1 , 1], since a(x) is there uniformly bounded away from zero. Hence, the solution propagator is ε-uniformly bounded (above and below) on [x 1 , 1], see Step 3 in the above proof. On the other hand, at the turning point x = 0 the Airy function has the value Ai(0), i.e. constant w.r.t. ε. Hence, the ε −1/6 scaling yields the asymptotic blow-up at x = 0.
This asymptotic blow-up is, for the time being, one of the key problems for extending asymptotic preserving schemes (like the WKB-method from [ABN11] or the adiabatic integrators from [LJL05]) up to the turning point. To mitigate this problem, yet still include the turning point into the scattering system, we made the simplifying assumption that a(x) ≡ x on some (small) interval [0, x 1 ]. This way we shall match the analytic solution on [0, x 1 ] to a numerical solution on [x 1 , 1]. The former is explicit up to a scaling factor that, however, inherits a numerical error from the approximation to ψ(1).

Numerical method and error analysis.
In this section we will review the WKB-marching method for the IVP (3.3) and derive error estimates for the BVP (2.4), extended to [0, 1].
We recall that the BVP (2.4) is solved in two steps: First the corresponding IVP (3.3) is solved numerically; then the numerical solution is scaled according to (3.4), to fit the right boundary condition. For the turning point problem we actually want to solve (2.4) on [0, 1]. But the solution to this BVP on [0, x 1 ] is given by a scaled Airy function as in (2.1). We are left with numerically approximating a solution on [x 1 , 1] and matching the two parts at x 1 to obtain a solution on the whole interval. In the 2-step solution process we incur an error from the WKB-marching method on [x 1 , 1] and this propagates into a second error from the (inaccurate) α-scaling of the "Airy-solution" on [0, The essential novelty compared to [AN18] is the inclusion of a first order turning point at x = 0. We proved in Section 4 that the exact solution ψ ε blows up at the turning point like Θ(ε − 1 6 ). Hence, one might expect that the corresponding numerical error would also be unbounded there. We recall that the numerical error stems from the α-scaling, where α(ψ,ψ ) depends on the numerically obtained approximations of ψ(1) andψ (1). In fact, the inclusion of the turning point "costs" a factor ε − 1 6 in the error estimates due to the blow-up of the sequence ψ ε at x = 0 (cp. the estimate (5.5) to (5.12) below). In these estimates the negative power of ε can be compensated by restricting the step size h in dependence of ε. With a turning point, this h(ε)-relation gets slightly less favorable.
In the special case of an explicitly integrable phase φ(x) from (2.7), the WKBmarching method is an asymptotically correct 3 scheme w.r.t. ε (with error order O(ε 3 )). This asymptotic correctness will compensate for the fact that the solution sequence ψ ε is unbounded at the turning point x = 0, and it will still yield an overall asymptotically correct scheme 4 for (2.4).
To clarify the notation, we summarize it in the following table. Here the superscript ( ) means that we refer to the function as well as its derivative. Let x 1 < x 2 < . . . < x N = 1 be a grid for the numerical method on [x 1 , 1], where x 1 > 0, and h := max 2≤n≤N |x n − x n−1 | is the step size. Occasionally we shall use a sub-or superscript ε to emphasize the ε-dependence of that quantity.
The Airy-WKB scheme for an approximation ψ Step 1: a) On [0, x 1 ]: The solution (2.1) satisfying the ICs in (3.3) reads 3 I.e. the numerical error decreases to zero with ε → 0, even for a fixed spatial grid. 4 Note: This constitutes a numerical scheme only for a coefficient function a(x) that is linear (or even quadratic, as in Section 6 below) on [0, x 1 ] for some x 1 ∈ (0, 1). b) On [x 1 , 1]: Compute a numerical approximation (ψ h,n , εψ h,n ) to the IVP (3.3) on the grid {x 1 , . . . , x N } via the WKB-marching method.
Step 2: We shall now review the basics of the second order WKB-marching method from [ABN11]: Here the focus is on the algorithm and error estimates. The background, including motivation for the tools used in this method can be read in [ABN11]. The method consists of two parts, first a transformation of the highly oscillatory problem (1.1) to a smoother problem, and second the numerical discretization of said smooth problem as to obtain an εasymptotically correct scheme.
Analytic transformation: The first order system (3.8) forŴ (x) is transformed to the system in the variable Z(x) as follows: where φ(x) is the phase function defined in (2.7). For this analytic transformation we assume that φ is explicitly available (for a generalization on numerically computed phases see Remark 5.3 below). This yields the system The above system exhibits much smoother solutions compared to the system forŴ (x) from (3.8). Moreover, the strong limit of its solutions Z ε as ε → 0 satisfies the trivial equation a numerical scheme that is at the same time ε-asymptotically correct and second order in the step size h. Numerical scheme: The second order (in h) scheme is rather non-standard and developed via the second order Picard approximation of (5.1): These (iterated) oscillatory integrals (with φ assumed to be known exactly) are then approximated using similar techniques as the asymptotic method in [INO06]. This yields the following scheme that is ε-asymptotically correct: For a given initial condition Z 1 := Z(x 1 ) the algorithm reads (5.2) Z n+1 = (I + A 1 n + A 2 n )Z n , n = 1, . . . , N − 1 , with the matrices A 1 n and A 2 n given as Here we used the following notations: In the end we obtain a sequence of vectors Z n which we have to transform back via This yields an approximation of the solution to the vector valued system (3.8) forŴ . Now let us formulate a discrete analogue of Lemma 3.1.
Proof. A proof can be carried out exactly as in [AN18, Lemma 3.5], with the only difference that the initial conditionŴ 1 is now ε-dependent (but uniformly bounded above and below).

Error estimates including the turning point.
The following result is a simple consequence of Theorem 3.1 in [ABN11]. It shows that the WKBmarching method applied to the IVP (3.3) -however, with a different IC compared to [ABN11]-yields the same h-and ε-order as when applied to the IVP proposed in [ABN11]. So we obtain: with a constant C independent of n, h and ε. Here, γ > 0 is the order of the chosen numerical integration method for evaluating the phase integral φ from (2.7).
Remark 5.3. In the error analysis of [ABN11], as well as in Proposition 5.2 above, the possible error of the matrices A 1 n , A 2 n in (5.2) arising from an incorrect phase φ is not taken into account. Errors of φ are only considered for the back transformation (5.3). However, in the recent, more complete error analysis in [AKU19], both occurrences of the error of φ are included. It would lead to a third error term in (5.5) that is O(ε 2 ) and includes the W 2,∞ -error of the phase (see [AKU19, Theorem 3.2]). We omit this additional error term here for brevity of the presentation.
Proof of Proposition 5.2. The proof is using the main result [ABN11, Theorem 3.1] for the IVP and it only remains to generalize it to the initial condition in (3.3). Let Y (x) = (y 1 (x), y 2 (x)) be the vector valued solution to (5.6) after the transformation via (3.6), and Y n is its numerical approximation at x n obtained via the WKB-marching method. Due to Theorem 3.1 of [ABN11] it holds Next we give a transformation formula to connect ϕ withψ, the solution of (3.3). One easily verifies that Using the asymptotic representations (3.2) we verify with C > 0 independent of n, h and ε.
In some applications the phase φ is exactly computable: In quantum tunneling models, e.g., the crystalline heterostructure leads to a piecewise linear potential V (x) (with jumps due to the contact potential difference), and hence piecewise linear a(x). In this case the h γ ε -error term in (5.5) drops out and the scheme satisfies an ε-uniform, second order in h error estimate. Moreover it is even asymptotically correct with respect to ε. The opposite situation, when φ(x) has to be computed numerically, will be discussed in Remark 5.5 below.
Scaling to fit the right boundary condition The numerical approximation W n for n = 1, . . . , N to the solution vector W (x n ) of (2.4) is obtained by first calculating the numerical approximationŴ n for the IVP (3.8) via the WKB-marching method. Then it is scaled withα :=α(Ŵ N ), i.e. W n :=αŴ n , n = 1, . . . , N , whereŴ N is the approximation toŴ (1) obtained in the last step of the WKB-scheme. Now we can give the error estimates for numerically solving the BVP (2.4) using thẽ α-scaled Airy functionα ψ − (with the choice c 0 := ε − 1 6 ) on [0, x 1 ] and theα-scaled numerical solutionαψ h,n on {x 1 , . . . , x N }, i.e. We recall that the ε-scaling of the Airy function on [0, x 1 ] is important here to satisfy the IC at x 1 for the IVP (3.3). Next we give error estimates for the hybrid solution (5.8), (5.9), i.e. Airy function on [0, x 1 ] coupled to the WKB-solution on [x 1 , 1]. While our main strategy follows §3.5 of [AN18], the turning point at x = 0, and thus, the unboundedness of ψ ε (0), causes technical challenges.
Theorem 5.4 (Convergence of the Airy-WKB method). Let Assumptions 2.1 and 2.3 be satisfied and 0 < ε ≤ ε 1 . Then the pair (ψ h , εψ h ) satisfies the following error estimates: a) In the region [0, x 1 ], we have |e h,n | + ε|e h,n | ≤ C h γ ε + Cε 3 h 2 , n = 1, . . . , N , where e h,n := ψ(x n ) −αψ h,n and e h,n := ψ (x n ) −αψ h,n . c) For the hybrid method on the interval [0, 1] we have the error estimate Remark 5.5. The h γ ε (or h γ ε 7/6 ) term drops out if the phase φ in (2.7) is explicitly integrable, leading to a second order scheme (in h) that is asymptotically correct with respect to ε. But when φ has to be computed numerically, e.g., via Simpson's rule where γ = 4, the scheme is still second order in h as long as h is bounded by O(ε is well known for the WKB approximation of highly oscillatory problems without a turning point (see [LJL05,ABN11]), yielding a scheme of order h 2 . Using a spectral method for the phase integral allows to drastically reduce the quadrature error for φ, as illustrated in [AKU19]. Proof.
[of Theorem 5.4] Within this proof, we will use Lemma 3.2 multiple times for argumentsŴ (x n ) as well asŴ n . Thus we choose δ := min(C 3 , C 5 ) with the lower bounds C 3 and C 5 on the arguments obtained in Lemma 3.1 and Lemma 5.1. It is crucial here that the solution vectorŴ (x n ) of the IVP (3.8) as well as the numerical approximationŴ n are in R 2 (see Lemma 5.1). Then the mapα is Lipschitz continuous with a constant Lα and uniformly bounded by a constant Cα, both independent of 0 < ε ≤ ε 0 , as stated in Lemma 3.2. a) For x ∈ [0, x 1 ] we first recall that cf. (3.4), (3.12). Hence, with (5.8) we have where we used the ε-uniform boundedness of the (scaled) Airy function Ai − x ε 2/3 on R + . We also note that the term ε − 1 6 cannot be compensated by Ai − x ε 2/3 , since the latter term takes a constant value (independent of ε) at x = 0. We also used the Lipschitz continuity ofα, in addition to Proposition 5.2. For the derivative we estimate as follows: This can be argued in the same manner as in Step 5 of the proof of Proposition 4.2. b) It is convenient to use the vector notation W on [x 1 , 1]. Note that scalingŴ with the constantα is equivalent to scaling (ψ, εψ ). Therefore the estimates after theα-scaling are In the second to last line we used the Lipschitz continuity and boundedness ofα as well as (3.11). In the last line we used the estimate (5.5) twice. Due to the (ε-uniform) equivalence of W and (ψ, εψ ) , the estimate above yields the desired bound on the interval [x 1 , 1]. c) The overall estimate on [0, 1] is a combination of the previous two.

Numerical results.
In this subsection we will present numerical results to illustrate the error estimates of Theorem 5.4 for the Airy-WKB method. Fig. 4 shows the error of the Airy-WKB method on [0, 1] with coefficient function a(x) = x and "switching-point" x 1 = 0.1. This is a convenient test case since its solution is explicitly known (an Airy function), yet the numerical WKB method on [x 1 , 1] is not trivial. Moreover, the integral of the phase φ(x), needed for the WKB-marching method, is explicitly available without numerical integration. For comparison we shall present two simulations, one with the exact phase φ(x) and one with a numerically computed phase (used both in (5.2) and (5.3)). In the right plot, the h γ ε 7/6 -term drops out of the error estimate (5.12), yielding an asymptotically correct scheme as ε → 0 (for fixed step size h).
In the left plot, the first term of (5.12) (i.e. h γ ε 7/6 , originating from the numerical integration of the phase φ) is dominant for large values of h and/or small values of ε. Hence, the error behaves like h 4 ε 7/6 , due to the Simpson rule with γ = 4, as visualized by the right slope triangle. At h = 0.1 we can clearly see an inversion of the error curves (ε = 2 −12 at the top and ε = 2 −8 at the bottom) with an ε-dependence of almost exactly O(ε −7/6 ). This error term originating from the numerical integration of the phase could be reduced to machine precision by using the spectral method proposed in [AKU19].
In the right plot (and likewise in the left plot for small h and large ε) we clearly observe the quadratic convergence rate in h for each value of ε. The error in ε is decreasing with order of about ε 3.4 to ε 3.6 and therefore better than the predicted estimates of order ε 17/6 from the second term in (5.12). This improved ε-order of the error does not originate in the choice of a linear potential, as a similar observation was already made in the error plot for the WKB-marching method in [ABN11, Fig.  3.1 (right)]. There, a simulation with a quadratic coefficient function was chosen and the error order in ε was showing better results than the predicted order of ε 3 .
For small values of both h and ε, the error is very small, such that it gets eventually polluted by round-off errors (due to double precision computations in Matlab). Hence, in this case the error is not showing a simple dependence on h and/or ε, and is around the order of 10 −11 .
In order to illustrate the efficiency of the Airy-WKB method, we make a comparison to a standard Runge-Kutta method, the 'ode45' single step Matlab solver. This solver is based on an explicit Runge-Kutta (4, 5) formula (the Dormand-Prince pair). It is adaptive to attempt to optimize the step size h using error estimates obtained from the comparison of a 4 th and 5 th order approximation. In this example we used (again) the linear coefficient function a(x) = x, such that we have the explicit formula for the exact solution at hand, in terms of the Airy function Ai. Additionally, the phase is explicitly available without numerical integration, and thus it does not contribute to an additional error or increase of the run time. Also, in this case the WKB scheme is asymptotically correct.
The goal of this test is to compare (for several fixed values of ε) the run times of the Airy-WKB and the Runge-Kutta methods, when both methods achieve (almost) the same numerical accuracy. To match the accuracies, we first determined the error of the Airy-WKB scheme, in the norm e h ∞ + ε e h ∞ (as defined in Theorem 5.4) for a grid with uniform step size h = 10 −3 . Then we specify an error tolerance for 'ode45' (somewhat by trial and error) to obtain a matching approximation error.  8.5833e-01 1.5724e-01 1.8105e-08 1.8443e-08 2 −10 3.2963e+00 1.6063e-01 1.3234e-09 1.2300e-09 run times stay constant for decreasing ε, since we leave the step size h fixed. But the Runge-Kutta method's run times grow strongly for smaller ε. This is a consequence of the following two contributing parts. One is due to the prescribed error tolerances. These need to go down, since the Airy-WKB method has decreasing error as ε → 0 (even for a fixed step size), and we need to obtain comparable errors for both methods. The second reason is that, for smaller ε, the oscillatory solution to the problem exhibits higher and higher frequencies, i.e. of order O( √ a ε ). Therefore the Runge-Kutta method needs to refine the step size h even more to resolve the oscillations and to meet the prescribed error tolerances.
6. Generalization to quadratic potentials close to the turning point. The goal of this section is a generalization of the hybrid method of Section 2-5 to the situation when the potential V (x) is quadratic instead of linear in the vicinity of the turning point (which is still of first order at x = 0). As we shall take a similar path as in the previous sections, not all of the (analogous) motivating deductions will be repeated. We start with specifying the assumptions on the coefficient function a(x) analogously to Assumption 2.1.
Assumption 6.1. Let parts a), b) and d) of Assumption 2.1 remain unchanged. c') More general than in Section 2-5, we now assume the potential to be quadratic in the left exterior and also in a (small) neighborhood of x = 0. More precisely we assume that ∃ x 1 ∈ (0, 1) such that a(x) = k 1 x 2 + k 2 x for x ≤ x 1 with k 1 < 0 and k 2 > −k 1 x 1 > 0 such that the second zero of a(x) is strictly larger than x 1 , and thus, not included in [0, x 1 ].
x V (x) E 0 expl. sol. x 1 numerical solution 1 expl. sol. Note that the above stated quadratic form of a(x) is in its most general form as a turning point at x = 0 requires a(0) = 0. We consider a potential V (x) → +∞ as x → −∞ (see Fig. 6) which requires a scattering solution for the equation ε 2 ψ +a(x)ψ = 0 to decay for x → −∞. For a quadratic potential, this solution is given on (−∞, x 1 ] by the parabolic cylinder function (PCF) denoted by U (ν, z), cf. [NHM10,§12].
For a(x) = k 1 x 2 + k 2 x, we have and a fundamental set of solutions is {U (ν, z(x)), U (−ν, iz(x))}. Here, U (ν, z(x)) is the solution that stays bounded (and even decays) for x → −∞. This yields the following BVP with transparent BCs: Here and in the sequel the notation for the derivative of the parabolic cylinder function is U (ν, z) where the always refers to the derivative w.r.t. the (second) argument z, not the order ν. An analogous result as Proposition 2.2 about existence and uniqueness of solutions holds for the BVP (6.2). As in Section 3 this BVP will first be reformulated as the IVP with a constant c 1 ∈ R \ {0} that we shall fix later. From here on we denote byψ the solution to this IVP. The solution ψ of (6.2) is then obtained by scalingψ to fit the right BC of (6.2), at x = 1. Note that this scaling preserves the validity of the left BC of (6.2), at x 1 .
6.1. Asymptotic blow-up at the turning point. In this section, we shall analyze the asymptotic behavior of solutions to the BVP (6.2) with quadratic potential close to the turning point. The following proposition will be used later to prove that the solution to (6.2) is not uniformly bounded w.r.t. ε ; this unboundedness arises at the turning point.
The lengthy proof is deferred to the Appendix A.3.
An essential requirement for the proofs in Section 3-5 was the uniform boundedness of the initial condition (ψ(x 1 ), εψ (x 1 )) w.r.t. ε. Proposition 6.2 shows that the initial condition in (6.3) has the asymptotic order O(g(µ)) w.r.t. ε, if c 1 is chosen independent of ε. In order to obtain again ε-uniform boundedness of the IC, we shall now choose c 1 depending on ε. As g(µ) is defined as an asymptotic series and g(µ) = Θ ε (h(µ)), we choose in the IC of (6.3): with h(µ) defined in (A.8). Using 1 h(µ) is more practical than 1 g(µ) since it can be implemented explicitly and it yields the same asymptotic scaling as 1 g(µ) . Then the IVP reads (6.5) After transforming (ψ, εψ ) toŴ ∈ R 2 via the matrix A(x) from (3.7) we get the vector-valued system (6.6) with the two matrices A 0 (x) and A 1 (x) as in (3.9).
Next we illustrate that the solutions of the BVP (6.2) become unbounded at the turning point x = 0, analogously to Example 4.1. First we consider Example 6.4. Consider (6.2) with x 1 = 0 and a(x) = x − x 2 2 for x ∈ [0, 1] and 0 < ε < 1. Then the explicit solution reads The quadratic coefficient a(x) = x − x 2 2 has a turning point at x = 0 and in Fig. 7 the absolute value of solutions |ψ ε (x)| is plotted for various ε > 0: The solutions |ψ ε (x)| are unbounded as ε → 0 at the turning point x = 0. Fig. 8 shows the plot of the family {ε|ψ ε (x)|}, which is bounded on [0, 1] uniformly in ε. It is even decreasing (in ε → 0) at x = 0.
The generalization of Example 6.4 to potentials as in Assumption 6.1 is the main result of the following proposition.
The lengthy and involved proof is provided in Appendix 11. 6.2. Numerical method and error analysis. In this section we extend the Airy-WKB method from Section 5 to the more general case of a quadratic potential in the vicinity of the turning point. As the fundamental solution to the equation in (6.2) for quadratic a(x) is a parabolic cylinder function (PCF), we will denote this method as PCF-WKB method. The error estimates and convergence results will be essentially the same as for the Airy-WKB method in Section 5.
Next we want to carry over the error estimates of the WKB-marching method [ABN11], just like in Proposition 5.2. Observe that, similar to (5.7), one can write the solutionψ(x) to the IVP (6.5) aŝ where ϕ is the solution to (5.6) and Using Proposition 6.2 one can see that [cos(ξ(t 1 )) − i sin(ξ(t 1 )) + O(ε)] = Θ ε (1) , and we get an analog error estimate as in Proposition 5.2: Under Assumptions 2.3 and 6.1, and for 0 < ε ≤ ε 1 it holds Here,Ŵ (x n ) is the exact solution to the IVP (6.6) at x n , andŴ n is its numerical approximation obtained by the WKB-marching method.
With this notation we can now formulate the analog result to Theorem 5.4 for quadratic potentials satisfying Assumption 6.1. The error orders are the same as in the case of a linear potential, since we considered a first order turning point at x = 0 in both cases.
Theorem 6.6 (Convergence of the PCF-WKB method). Let Assumptions 2.3 and 6.1 be satisfied and 0 < ε ≤ ε 1 . Then the pair (ψ h , εψ h ) satisfies the same error estimates as in Theorem 5.4.
The h γ ε 7/6 term from the numerical integration of the phase φ(x) from (2.7) appears again in the error estimates of Theorem 6.6. As mentioned in Section 5.2, this term drops out for some applications, when the phase is explicitly integrable.
6.3. Numerical results. The illustration of the error estimates of Theorem 6.6 for the PCF-WKB method is done analogously to Subsection 5.3. Here the coefficient function of (1.1) is chosen as a(x) = x − x 2 2 on [0, 1]. Therefore the solution ψ is explicitly known as a parabolic cylinder function. The calculations are done in MATLAB R , where the PCF is not implemented, but can be obtained via relations to other functions, i.e. the confluent hypergeometric function [NHM10, §12.7(iv)]. As the range of ε we chose ε = 2 −5 , . . . , 2 −9 , since for ε = 2 −10 the evaluation of the PCF returns Inf as values. This is because its order ν becomes large (negative) for small ε, and thus, the PCF maps to very large values. For ε = 2 −9 the evaluation already gets very inaccurate, such that we used Mathematica R for this case -to evaluate the PCF for the reference solution.
For the plots in Fig. 9, we computed the integral for the phase φ(x) numerically, once for the plot on the right with an error tolerance of 10 −12 (such that the corresponding error term is negligible in comparison to the second error term in (5.12)), and once using the composite Simpson rule for the phase in the plot on the left. We clearly see the h 2 convergence rate of the method for each ε. An interesting difference can be seen in the left plot for small ε and large h: Here the h 4 ε 7/6 -term, originating from the numerical integration of the phase, is dominant.
In the region with quadratic convergence (i.e. when the second term in (5.12) is dominant), the error in ε is decreasing with order of about ε 3 to ε 3.7 , and thus it is again (cf. Subsection 5.3) slightly superior to the predicted estimates of order ε 17 6 .
7. Generalization and outlook. The next natural question is how to alleviate the restriction that the potential V (x) should be (exactly) linear or quadratic in a neighborhood of the turning point. An obvious and frequently used strategy is to approximate the potential: In [HH08] a piecewise constant approximation was  But, as we shall illustrate next, such linear (or even higher order) approximations lead to errors that are unbounded as ε → 0. Hence, they cannot serve as a starting point to construct uniformly accurate schemes. The error encountered by taking Airy functions (as solutions to the reduced problem -i.e. the Airy equation) can only be uniformly bounded when confining to a region around the turning point that shrinks fast enough as ε → 0 (particularly like o(ε 2/5 ), [Mil06]). But, for the time being, the WKB-marching method requires a constant-in-x transition point x 1 . Next we illustrate the consequence of a linear approximation of the potential close to the turning point: We chose a piecewise quadratic coefficient function a(x), see Fig.  10, along with its linear approximationã on [0, x 1 ] with some (fixed) x 1 ∈ (0, 1).
Since the solutions of the scattering problem for both potentials are analytically known in terms of Airy functions and parabolic cylinder functions, the error plotted in Fig. 11 is not due to any numerical method. The non-uniform error (in ε) stems only from the modified coefficient function, and hence modified phase of the solution: For ε small, the approximate solution is entirely out of phase close to x = x 1 . On the one hand the error in Fig. 11 is of order O(x 3 1 ) for fixed ε, which stems from the linear approximation of the potential. On the other hand, for fixed x 1 the error grows approximately like ε −3/2 . As a conclusion of this feasibility study and as an outlook for a follow-up work there are two options to proceed. Since a linear approximation of a(x) leads to εuniform errors for transition points x 1 = o(ε 2/5 ) (see [Mil06]), a first option would be to extend the WKB-marching method to ε-dependent intervals of the form [ε α , 1], with some 0 < α < 1. Since the WKB-method yields an O(ε 3 )-error (for analytically integrable phase functions, see Proposition 5.2) on ε-independent intervals, an exten-  x 3 1 Fig. 11: Error of the scattering solution due to the linear approximation of a(x) byã(x), as shown in Fig. 10. The plot shows the error in dependence of ε and the cut-off point x 1 ; apparently it is of the order x 3 1 /ε 1.5 . The function E(x) is the difference of the two solutions.
A second option is to use Langer functions [Lan31] instead of Airy functions on a fixed interval [0, x 1 ], coupled to the WKB-marching method. Langer functions are essentially Airy functions with a (generally) non-linear transformation of the argument and a (linear) modification of the amplitude, cf. modified Airy functions in The explicit definitions of A s (ζ), B s (ζ), C s (ζ) and D s (ζ) are not needed here and can be found at [NHM10,12.10.42(44)]. We just note that they are independent of ε and that (A.10) A 0 (ζ) = 1 ; D 0 (ζ) = 1 .
These asymptotic expansions are in terms of the Airy function and their validity region includes one (x = 0 ⇔ t = 1) of the two turning points for the quadratic coefficient function. On top of that they hold uniformly for all t ∈ [−1 +δ, 1], i.e. bounded away from the second turning point at t = −1 by some arbitrarily small δ > 0. Within the region [−1 +δ, 1] the coefficient function 1 − t 2 of the transformed equation is non-negative. In order to keep the notation clear, we shall mostly stick to the notation from the literature, i.e. using µ and ζ(t).
We also need some properties for the function ϕ(t): asymptotic expansion. But if two asymptotic expansions are power series (so called Poincaré-type expansions), their product is again an asymptotic expansion (see [Ol74, Ch.1, §8.1(ii)]).
The asymptotic expansions we are interested in are power series with respect to the asymptotic sequence ε 2k . For some f ε (t) ∼ ∞ k=0 a k (t)ε 2k and g ε (t) ∼ it holds that (f ε · g ε )(t) ∼ We shall now apply this to the product Ai(µ Here ζ(t) is an ε-independent function in t, but η(t) is not only t-dependent but also of order O(ε −1 ) as ε → 0. As η(t) appears only within sin and cos, it does not affect the ε-asymptotic behavior of the representation.
In the same manner we get With Proposition 6.2 we can now prove Proposition 6.5.
Proof of Proposition 6.5. For readability of this proof, we omit the index ε in ψ ε (the solution to the BVP (6.2)) andψ ε (the solution to the IVP (6.5)). The proof of Proposition 4.2 relies on the fact that the scaled Airy function solution on [0, x 1 ] only depends on the variable −xε −2/3 (see Remark 4.3). Since this is not true for the PCF solution, the strategy of this proof will deviate from Proposition 4.2, and we shall also need the asymptotic expansion of U (ν, z(x)). Still, we make a case distinction similar to the proof of Proposition 4.2: Region x 1 ≤ x ≤ 1: In this region the function a(x) is C 2 and satisfies a(x) ≥ τ 1 > 0. Hence, Lemma 6.3 yields Ŵ (x) = Θ ε (1) on [x 1 , 1] for the vector valued solutionŴ to the IVP (6.6). Since the norms ofŴ (x) and (ψ(x), εψ (x)) are (ε-uniformly) equivalent, we get ψ (x) εψ (x) = Θ ε (1) , x 1 ≤ x ≤ 1 .
For the vector solution of the BVP (6.2) in the region x 1 ≤ x ≤ 1 this yields whereψ(x) solves the IVP (6.5).
We start with the proof of statement a). For the parabolic cylinder function U (ν, z(x)) we can use the asymptotic expansion (A.6) to get the following asymptotic representation (when using only the first term with A 0 (ζ) = 1): (A.21) U (ν, z(x)) = g(µ) 2 √ πϕ(t)µ