3.1. Construction of the Stokes diagram

Stokes diagram provides a graphical representation illustrating how WKB solutions change dominance across different regions in the complex $\mathbb{C}$ plane and the detailed construction includes the following. (i) Determining the turning points $x_t$ by solving $Q(x_t)=0$ . (ii) Plotting the Stokes lines emanating from each turning point, defined by $\text{Im}(\int _{x_t}^x S_{-1}\,{\rm d}x)=0$ . Along a Stokes line, the WKB wave functions exhibit purely exponential behaviour without phase variation. The relative dominance between the two eigenfunctions is determined by the sign of $\text{Re}(\int _{x_t}^xS_{-1}\,{\rm d}x)$ . We use outward (inward) arrows to indicate the ‘+’ (‘–’) eigenmode dominance. (iii) Introducing branch cuts to partition the Riemann surface into sheets mapped bijectively to the $\mathbb{C}$ plane.

3.2. Exact WKB eigenfunctions

The solution of (2.3) can be expressed as a linear combination of exact WKB eigenfunctions $\phi ^\pm$ , i.e. $\phi (x)=c^+\phi ^+(x)+c^-\phi ^-(x)$ , where $c^\pm$ are constant coefficients.

In the vicinity of a simple regular turning point satisfying $\partial _xQ(x_t)\neq 0$ , the original ODE can be mapped to the canonical Airy equation $\delta ^2 \partial _z^2\psi (z)-z\psi (z)=0$ via a local coordinate transformation $z(x)$ expanded as $z(x)\equiv \sum _{n,m}\delta ^n z_{n,m}\Delta x^m$ , where $\Delta x \equiv x-x_t$ . The practical steps for computing the expansion coefficients $z_{n,m}$ are given by Kawai & Takei (Reference Kawai and Takei2005, Chapter 2.3). The exact WKB eigenfunctions in this region are

(3.4) \begin{align} \phi ^\pm =\psi ^\pm (\mathcal{L}^{-\text{1}}\mathcal{B}[z(x)]). \end{align}

Away from the turning point, the eigenfunctions take the form

(3.5) \begin{align} \phi ^{\pm }=\psi ^\pm \left (\mathcal{L}^{-\text{1}}\mathcal{B}[z(x_0)]\right )\exp \left (\int _{x_0}^x\mathcal{L}^{-\text{1}}\mathcal{B}[S^{\pm }]\,{\rm d}x\right )\!, \end{align}

where $x_0$ is a reference point chosen sufficiently far away from the turning point that both the Padé approximation of the Borel-resummed quantity $\mathcal{B}[S^{\pm }]$ and Taylor expansion of $z(x)$ in $\Delta x$ remain valid.

3.3. Connection formula

The solution path, i.e. the daughter-wave trajectory defined on the complex $\mathbb{C}$ plane, is typically taken along the real axis. We denote its initial and final points by $x_a$ and $x_b$ , respectively.

In contrast to the conventional analytical approach, coefficients $c^\pm$ are constant only within each Stokes sector. Due to the discontinuous behaviour of the exact-WKB eigenfunctions across the Stokes lines, $c^\pm$ must undergo corresponding discontinuous jumps to ensure the continuity of the full solution $\phi (x)$ .

The connection formula illustrates how the coefficients $c^\pm$ are transformed along the path. We denote by $c^\pm _a$ and $c^\pm _b$ the values of $c^\pm$ at the endpoints $x_a$ and $x_b$ , respectively. The connection formula is then expressed as

(3.6) \begin{align} \left (\begin{array}{c} c_{b}^{+}\\[3pt] c_{b}^{-} \end{array}\right )= \unicode{x1D648}_{j}\cdots \unicode{x1D648}_{2} \unicode{x1D648}_{1}\left (\begin{array}{c} c_{a}^{+}\\[3pt] c_a^- \end{array}\right ). \end{align}

The matrices $ \unicode{x1D648}_{1\sim j}$ represent the transition matrices induced each time the solution path crosses the Stokes lines and branch cuts, encoding the jump of $c^{\pm }$ . The explicit form of $ \unicode{x1D648}$ depends on the local crossing geometry between the solution path and the lines. In particular, when the solution path crosses a Stokes line, the transition matrix is $ \unicode{x1D648}_{\pm }=\begin{pmatrix} 1 & \quad i\\ 0 & \quad 1 \end{pmatrix}$ for Stokes lines where ‘+’ eigenmode is dominant and $\begin{pmatrix} 1 & \quad 0\\ i & \quad 1 \end{pmatrix}$ for Stokes lines where ‘–’ eigenmode is dominant. By default, the wave path is taken counter-clockwise around the relevant turning point; for clockwise traversal, the corresponding transition matrix is $ \unicode{x1D648}_{\pm }^{-1}$ . Similarly, crossing a branch cut corresponds to the transition matrix $ \unicode{x1D648}_c=\begin{pmatrix} 0 & \quad -i\\ -i & \quad 0 \end{pmatrix}$ .

When a wave sequentially encounters Stokes lines from different turning points, a normalisation matrix must be applied before introducing the transition matrix of the subsequent turning point:

(3.7) \begin{align} \unicode{x1D648}_V=\left (\begin{array}{c@{\quad}c} \mathcal{L}^{-\text{1}}\mathcal{B}[e^{({1}/{2})\oint _I S^{+}{\rm d}x}] & 0\\[5pt] 0 & \mathcal{L}^{-\text{1}}\mathcal{B}\left [e^{({1}/{2})\oint _{I}S^{-}{\rm d}x}\right ] \end{array}\right ). \end{align}

The total transition matrix is then given by $ \unicode{x1D648}_V^{-1} \unicode{x1D648} \unicode{x1D648}_V$ , with $ \unicode{x1D648}$ in $ \unicode{x1D648}_\pm$ , $ \unicode{x1D648}_c$ . This matrix ensures the consistency of the global connection procedure. The Voros multiplier $V\equiv ({1}/{2})\oint _I(S^+-S^-)\,{\rm d}x$ , with $I$ denoting the closed contour encircling two associated turning points, is a key quantity for analysing wave coupling and potential well/barrier problems. Methods for calculating the Voros multiplier can be found elsewhere (Kreshchuk & Gulden Reference Kreshchuk and Gulden2019; Grassi, Hao & Neitzke Reference Grassi, Hao and Neitzke2022).

3.4. Boundary condition of a daughter wave

For convective PIs, one is often interested in how a pure incident daughter-wave perturbation at $x=x_a$ nonlinearly grows to a significant amplitude when it reaches a probe point $x=x_b$ . In conventional WKB, one might set initial conditions such as $\phi (x_a)=c_a, \phi ^\prime (x_a)=c_a S^+_{-1}(x_a)/\delta$ to represent a pure incident ‘+’ eigenmode. However, this introduces WKB approximation errors, generating a small admixture of the ‘–’ eigenmode that might grow exponentially and obscure the ‘+’ eigenmode of interest.

In the exact WKB framework, a pure incident daughter wave can be represented without admixture by

(3.8) \begin{align} \phi (x_a)=c^+_a, \phi ^\prime (x_a)=c_a^+ \mathcal{L}^{-\text{1}}\mathcal{B}[S^+(x_a)]. \end{align}

During the nonlinear coupling, mode exchange can occur: a daughter wave corresponding to ‘+’ eigenmode at $x_a$ may become a ‘–’ eigenfunction at $x_b$ , or vice versa. Considering this scenario, the amplification factor of this daughter wave is

(3.9) \begin{align} A_{\text{inhom}}=\ln \dfrac {|c_b^-\phi ^-(x_b)|}{|c_a^+\phi ^+(x_a)|}=\ln \frac {|c_b^-|}{|c_a^+|} +\ln \dfrac {|\phi ^{-}(x_b)|}{|\phi ^{+}(x_a)|}, \end{align}

with the contributions from both the Stokes jump of constant coefficients and the amplitude difference of the eigenmodes at the endpoints, given by the two terms on the right-hand side of (3.9), respectively.

4.1. Resonant forward scattering

In the resonant forward scattering regime, both daughter waves ‘1’ and ‘u’ are resonant and propagate from $-\infty$ to $+\infty$ . As shown in figure 2(a), wave ‘1’ corresponds to eigenmode ‘–’ at $x\to -\infty$ and ‘+’ at ‘ $x\to +\infty$ ’, while wave ‘u’ behaves inversely. Assuming the coefficients $c_a^\pm$ before coupling are known, the coefficients $c_b^\mp$ after coupling can be determined by the connection formula.

According to the Stokes diagram in figure 1(a), the wave encounters two Stokes lines and the connection formula is

(4.5) \begin{align} \left (\begin{array}{c} c_{b}^{+}\\[3pt] c_{b}^{-} \end{array}\right )=\left (\begin{array}{c@{\quad}c} 1 & 0\\[3pt] ie^{2V/\epsilon } & 1 \end{array}\right )\left (\begin{array}{c@{\quad}c} 1 & -i\\[3pt] 0 & 1 \end{array}\right )\left (\begin{array}{c} c_{a}^{+}\\[3pt] c_{a}^{-} \end{array}\right ). \end{align}

Consequently, for wave ‘1’, the Stokes jump yields

(4.6) \begin{align} c_{a}^{-}\to c_{b}^{+}=c_{a}^{+}-ic_{a}^{-}. \end{align}

However, for wave ‘u’, the jump result yields

(4.7) \begin{align} c_{a}^{+}\to c_{b}^{-}=ic_{a}^{+}e^{2V/\epsilon }+c_{a}^{-}\left (1+e^{2V/\epsilon }\right ). \end{align}

The final amplitude of one daughter wave depends not only on its incident amplitude, but also on that of the other daughter wave. Since daughter waves are typically seeded by background noise, arbitrary initial conditions $c_a^\pm$ may in principle be realised. Special cases, such as $c_a^+=ic_a^-$ , which lead to the vanishing of wave ‘1’ after coupling, are also possible. However, our physical interest usually lies in the strongly amplified daughter waves that grow from the noise floor to amplitudes large enough to affect the system dynamics. Since this requires amplification over many orders of magnitude, the relevant regime is typically of large amplification ( $A\gg 1$ ). In this limit, even a modest difference in amplification factor is exponentially magnified, such that the seed corresponding to the maximum amplification overwhelmingly dominates the final wave amplitude. Therefore, for most practical purposes, such as evaluating the PI-induced lower hybrid sidebands in tokamaks (Cesario et al. Reference Cesario, Cardinali, Castaldo, Paoletti and Mazon2004) or estimating the worst gain induced by stimulated Raman scattering in inertial confinement fusion (Rosenberg et al. Reference Rosenberg2018), it is sufficient to consider only the boundary condition that yields the maximum amplification.

In figure 4, the actual wave amplitude evolution is traced by combining the eigenmode evolution with the Stokes jumps. When focusing on the evolution of a single wave, the incident coefficient of the other daughter wave is set to zero. It can be observed that, although the final amplification predicted by the exact WKB and the conventional WKB approaches is identical, the local wave behaviours exhibit clear differences near the boundaries of the coupling region, i.e. in the vicinity of the turning points. In this region, higher-order WKB corrections become significant and the conventional WKB validity condition breaks down. Remarkably, the higher-order contributions from the two turning points cancel each other, which ultimately leads to the same prediction for the final amplification factor.

Figure 4.

Amplitude evolution of (a) daughter wave ‘u’ and (b) daughter wave ‘1’ for resonant forward scattering PI. From the exact WKB perspective, the green solid lines represent waves with well-defined wavenumbers, while purple-shaded regions indicate intervals where the daughter waves are strongly coupled and cannot be clearly distinguished. The blue dashed lines show the amplitude evolution predicted by the conventional WKB method. The incident coefficients are $c_a^+=1, c_a^-=0$ for panel (a) and $c_a^+=0, c_a^-=1$ for panel (b).

Consistency checks of the exact WKB solutions of the original ODE (2.3) are presented in figure 5, where they are found to be in indistinguishable agreement with numerical solutions and the parabolic cylinder representation.

Figure 5.

The $\phi _1$ solution of (2.3) obtained via exact WKB (solid line), a numerical ODE solver (crosses) and the parabolic cylinder function representation in (4.2) (circles) for the resonant forward scattering regime. Black curves/markers denote the amplitude of $\phi _1$ , while red curves/markers denote the unwrapped phase. Parameters of panels (a) and (b) are the same as those used in figures 4(a) and 4(b), respectively.

4.2. Resonant backward scattering

In the resonant backward scattering regime, wave ‘u’ propagates in the opposite direction. From the wave-vector evolution shown in figure 2(b), we identify that the ‘+’ eigenmode at $x\to -\infty$ and the ‘–’ eigenmode $x\to +\infty$ correspond to wave ‘1’, whereas the ‘–’ eigenmode at $x\to -\infty$ and ‘+’ eigenmode at $x\to +\infty$ correspond to wave ‘u’. Therefore, the two coefficients representing the incident waves are $c_a^+$ and $c_b^+$ . Indicated by the Stokes diagram in figure 1(b), the connection formula is

(4.8) \begin{align} \left (\begin{array}{c} c_{b}^{+}\\[3pt] c_{b}^{-} \end{array}\right )=\left (\begin{array}{c@{\quad}c} 1 & ie^{2V/\epsilon }\\[3pt] 0 & 1 \end{array}\right )\left (\begin{array}{c@{\quad}c} 1 & 0\\[3pt] -i & 1 \end{array}\right )\left (\begin{array}{c} c_{a}^{+}\\[3pt] c_{a}^{-} \end{array}\right )\!. \end{align}

For wave ‘1’, the Stokes jump is

(4.9) \begin{align} c_{a}^{+}\to c_{b}^{-}=ie^{-2V/\epsilon }(c_{a}^{+}-c_{b}^{+}), \end{align}

while for wave ‘u’, which propagates in the opposite direction to wave ‘1’, the Stokes jump reads

(4.10) \begin{align} c_{b}^{+}\to c_{a}^{-}=ie^{-2V/\epsilon }(c_{a}^{+}-c_{b}^{+})+ic_{a}^{+}. \end{align}

Note that in backward-scattering regimes, $V_r\lt 0$ because $v_{gu}v_{g1}\lt 0$ . Since we are only interested in the maximally amplified PI response, we focus on cases in which the large factor $\exp (2|V_r|/\epsilon )$ appears in the post-coupling coefficients $c_a^-$ and $c_b^-$ . Combinations such as $c_a^+=c_b^+$ are therefore disregarded.

The amplitudes of the $\pm$ eigenmodes shown in figure 3(b) vary little across space, yet the amplitude difference between the ‘+’ and ‘–’ eigenmodes becomes pronounced. This disparity arises from the $\mathcal{L}^{-\text{1}}\mathcal{B}[S^+]$ integral evaluated from the upper turning point in figure 1(b) to 0, together with the $\mathcal{S}[R^-]$ integral from the lower turning point to 0. These two contribute to an amplitude disparity approximately $|\phi _-/\phi _+|\simeq \exp (-|V_r|/\epsilon )$ .

Taking into account that the Stokes jump already provides an amplification of $\exp (-2V_r/\epsilon )$ , the final amplification factors of both daughter waves are estimated as $\exp (-V_r/\epsilon )$ . To confirm this precisely, figure 6 shows the amplitude evolution from the perspective of exact WKB. Again, when examining the growth of a particular daughter wave, the incident coefficient of the other wave is set to 0. The logarithmic amplitude difference between ends of the interval results in $7.854$ for both daughter waves, consistent with the estimated amplification factor $A_{\text{inhom}}=|V_r|/\epsilon$ , identical to the resonant forward scattering case. Verification of the exact WKB solution for this regime is shown in figure 7, demonstrating full consistency among all methods.

Figure 6.

Amplitude evolution of (a) wave ‘1’ and (b) wave ‘u’ in resonant backward scattering PI. The purple shaded regions represent coupling intervals where the two daughter waves cannot be distinguished from each other. The incident coefficients are chosen as $c_a^+=1, c_b^+=0$ in panel (a) and $c_a^+=0, c_b^+=1$ in panel (b).

Figure 7.

Comparison of the solution $\phi _1$ of (2.3) obtained from the exact-WKB method (solid lines), a numerical ODE solver (crosses) and the parabolic-cylinder-function representation (4.2) (circles) for the resonant backward scattering regime. Black curves/markers denote the amplitude of $\phi _1$ , and red curves/markers denote the unwrapped phase. Panels (a) and (b) use the same parameters as figures 6(a) and 6(b), respectively.

The $\exp (|V_r|/\epsilon )$ amplification appears to contradict the zero amplification predicted by (2.5). From the conventional WKB perspective, the zero amplification results from the fact that none of the phase integrands $\pm \sqrt {Q(x)}$ on the entire real axis possesses a real part, except for small corrections of order $O(\epsilon )$ .

The discrepancy originates from an intrinsic limitation of the conventional WKB method, which ignores Stokes jumps. This ignorance of Stokes jumps is equivalent to an implicit boundary condition that the amplitude of the unconcerned wave vanishes at the incident side of the concerned wave; only then can the eigenmode corresponding to the incident concerned wave evolve without a coefficient jump. In the resonant forward scattering regime, where both daughter waves propagate in the same direction, this assumption causes no issue. However, in the resonant backward scattering regime, it effectively sets the outgoing amplitude of the counter-propagating wave ‘u’ to zero. Since parametric amplification relies on a positive-feedback loop between both waves, suppressing one prevents the growth of the other, leading to the zero amplification predicted by conventional WKB. This can be verified by substituting $c_a^-=0$ into (4.9) and (4.10), which yields $c_b^-=-ic_a^+$ and a final attenuation of wave ‘1’ about $\exp (-|V_r|/\epsilon )$ .

As shown within the exact WKB framework, the zero amplification predicted by the conventional WKB treatment for resonant backward scattering PIs ( $({\varGamma \varGamma ^*}/{v_{gu}v_{g1})}\lt 0$ ) is invalid, due to the omission of boundary-condition data, in particular, Stokes jumps. In our previous work (Zhou et al. Reference Zhou, Xiang, Gan and Xia2025), the condition $\alpha \alpha ^*\lt 0$ was expressed in terms of $\alpha =\varGamma /v_{gu}$ and $\alpha ^*=\varGamma ^*/v_{g1}$ without distinguishing whether it arises from $\varGamma \varGamma ^*\gt 0$ , $v_{gu}v_{g1}\lt 0$ or $\varGamma \varGamma ^*\lt 0$ , $v_{gu}v_{g1}\gt 0$ (in the latter case, $\varGamma ^*$ is no longer the complex conjugate of $\varGamma$ ). The present analysis shows that these cases lead to fundamentally different behaviours. The former, i.e. the resonant backward scattering considered here, yields finite amplification, whereas the latter corresponds to the resonant forward scattering with zero amplification. Nevertheless, the previously reported mechanism, whereby increasing the linear growth rate of the daughter wave ‘u’ in the latter regime breaks the negative feedback loop and restores PI amplification, remains applicable.

4.3. Non-resonant forward scattering

Compared with the resonant type, where the incident amplitudes of both daughter waves must be prescribed, the non-resonant regime allows a much simpler boundary condition, as the heavily damping quasi-mode cannot exist as an incident wave and is instead generated exclusively through nonlinear coupling.

Figure 2(c) shows that the ‘+’ eigenmode corresponds to daughter wave ‘1’ and the ‘–’ eigenmode corresponds to the quasi-mode ‘u’ throughout the domain, without any mode exchange as in the resonant regimes.

Given that the quasi-mode’s incident coefficient $c_a^-$ is zero and Stokes diagram in figure 1(c), we have the connection formula:

(4.11) \begin{align} \left (\begin{array}{c} c_{b}^{+}\\[3pt] c_{b}^{-} \end{array}\right )=\left (\begin{array}{c@{\quad}c} 1 & 0\\[3pt] ie^{2V/\epsilon } & 1 \end{array}\right )\left (\begin{array}{c@{\quad}c} 1 & 0\\[3pt] i & 1 \end{array}\right )\left (\begin{array}{c@{\quad}c} 1 & i\\[3pt] 0 & 1 \end{array}\right )\left (\begin{array}{c} c_{a}^{+}\\[3pt] 0 \end{array}\right ). \end{align}

Solving the connection formula, we find that the Stokes jump for wave ‘1’ is

(4.12) \begin{align} c_{a}^{+}\to c_{b}^{+}=c_{a}^{+}, \end{align}

while the quasi-mode’s Stokes jump yields

(4.13) \begin{align} c_{a}^{-}=0\to c_{b}^{-}=ic_{a}^{+}\left (1+e^{2V/\epsilon }\right ). \end{align}

Therefore, amplification of wave ‘1’ does not depend on the Stokes jump and is entirely governed by its intrinsic spatial evolution, as illustrated in figure 8(a). Because of the absence of Stokes jump contribution and mode exchange, the conventional WKB treatment remains valid in this regime.

Figure 8.

Amplitude evolution of (a) wave ‘1’ and (b) quasi-mode ‘u’ for non-resonant parametric instability (PI) within the exact WKB framework. The behaviour of wave ‘1’ in panel (a) is identical for forward and backward scattering; thus, only one case is shown. The conventional WKB result is included as a blue solid line for comparison. In panel (b), the quasi-mode ‘u’ exists only over part of the spatial interval, as it is generated solely through nonlinear coupling.

In the conventional WKB approach, increasing the linear spatial damping rate $K_{ui}$ of quasi-mode ‘u’ causes the effective growth rate to spread across the entire real axis, rather than remaining localised between the turning points. Consequently, the integral evaluated over a finite region becomes smaller than the standard value $V_r/\epsilon$ .

We also compare the conventional WKB prediction with the evolution obtained from the exact WKB analysis in figure 8(a). The final amplification factor of daughter wave ‘1’ between a fixed spatial interval is found to be identical in both approaches. Although the higher-order WKB terms included in the exact WKB framework modify the local amplitude, these contributions are axially symmetric about the midpoint between the two turning points, and thus their net contribution to the integral over the entire interval vanishes.

The evolution of the quasi-mode is shown in figure 8(b). It is triggered by a Stokes jump at the second crossing point, with an amplitude much smaller than wave ‘1’, and quickly decays due to strong linear damping. This may account for the difficulty in experimentally detecting the quasi-mode. Figure 9(a) compares exact WKB solutions in this regime with other methods, showing that the results of all the approaches are numerically distinguishable.

Figure 9.

Comparison of the solution $\phi _1$ of (2.3) obtained from the exact-WKB method (solid lines), a numerical ODE solver (crosses) and the parabolic-cylinder-function representation in (4.2) (circles). Panel (a) is the non-resonant forward scattering PI regime, while panel (b) is the non-resonant backward scattering. Black curves/markers denote the amplitude of $\phi _1$ and red curves/markers denote the unwrapped phase. The linear growth trend observed in the parabolic cylinder representation and numerical solution in panel (b) indicates the onset of numerical stability.

4.4. Non-resonant backward scattering

Finally, the results in the non-resonant backward scattering case are nearly identical to those of the non-resonant forward scattering case. The key difference is that the quasi-mode ‘u’ propagates and damps in the negative $x$ -direction, leading to the boundary condition $c_b^-=0$ .

The connection formula based on Stokes diagram in figure 1(d) is

(4.14) \begin{align} \left (\begin{array}{c} c_{b}^{+}\\[3pt] 0 \end{array}\right )=\left (\begin{array}{c@{\quad}c} 1 & -ie^{2V/\epsilon }\\[3pt] 0 & 1 \end{array}\right )\left (\begin{array}{c@{\quad}c} 1 & 0\\[3pt] -ie^{-2V/\epsilon } & 1 \end{array}\right )\left (\begin{array}{c@{\quad}c} 1 & 0\\[3pt] -i & 1 \end{array}\right )\left (\begin{array}{c} c_{a}^{+}\\[3pt] c_{a}^{-} \end{array}\right ). \end{align}

The Stokes jump of wave ‘1’ is the same as (4.12), indicating no contribution to PI amplification. The Stokes jump of quasi-mode ‘u’ is given by

(4.15) \begin{align} c_b^-=0\to c_a^-=ic_a^+(1+e^{-2V/\epsilon }). \end{align}

The intrinsic amplitude evolution of wave ‘1’ is the same as in the non-resonant forward regime shown in figure 8(a). As shown in figure 8(b), the quasi-mode ‘u’ is generated at the crossing point $P_2$ in figure 1(d), with again an amplitude that remains negligibly small and rapidly decays due to linear damping.

Cross-validation for this regime is presented in figure 9(b). The numerical solution becomes unstable beyond $\text{Re}(x)=-0.74$ , as the quasi-mode exponentially amplifies the small numerical errors during propagation toward increasing x. The analytical solution based on parabolic cylinder functions also breaks down beyond $\text{Re}(x)=1.08$ , where severe cancellation occurs between the two large contributions satisfying $c^{[1]}\phi _1^{[1]}\simeq -c^{[2]}\phi ^{[2]}_1$ . In contrast, the exact WKB solution remains stable throughout the entire domain, maintaining a well-behaved amplitude and phase without any sign of breakdown. This demonstrates the superior robustness of the exact WKB approach compared with both direct numerical integration and analytical representations via special functions.