The WentzelKramersBrillouin (WKB) approximation is frequently used to explore the mechanics of the cochlea. As opposed to numerical strategies, the WKB approximation facilitates analysis of model results through interpretable closedform equations and can be implemented with relative ease. As a result, it has maintained relevance in the study of cochlear mechanics for half of a century. Over this time, it has been employed to study a variety of phenomena, including the limits of frequency tuning, active displacement amplification within the organ of Corti, feedforward mechanisms in the cochlea, and otoacoustic emissions. Despite this ubiquity, it is challenging to find rigorous exposition of the WKB approximation's formulation, derivation, and implementation in cochlear mechanics literature. In this tutorial, the foundations of the WKB approximation are discussed in application to models of one and twodimensional cochlear macromechanics. This includes mathematical background, rigorous derivation and details of its implementation in software.
I. INTRODUCTION
Models of one (1D) and twodimensional (2D) macromechanics have offered some of the most significant interpretations of cochlear physiology, historic and modern. It is intuitive that threedimensional (3D) models should offer more physically realistic results than 2D or 1D models, but this alone implies a “more the merrier” view of model dimensionality that coincides with quantitative accuracy but not with frequency of implementation or impact on the field of cochlear mechanics.
Important results of 1D models include the existence and character of stapesdriven traveling waves and the presence of a region of negative damping,^{1–7} as well as intracochlear reflections and otoacoustic emissions (OAEs).^{5,8–15} Qualitative similarity across frequency/space and quantitative similarity in the longwave region to in vivo cochlear responses make 1D models attractive for the exploration of fundamental macromechanical phenomena. Implementation and modification of the dynamics to account for features such as nonlinearity and roughness are also far simpler in 1D models than in 2D or 3D models.^{10–13}
On the other hand, 2D models allow for more physical results in the shortwave and cutoff regions of the cochlear response than 1D models,^{16–21} allowing for more complete exploration of potential modecoupling phenomena in the traveling wave.^{22–24} Moreover, 2D models are able to capture fluid mechanical properties in the scalae, which allows for interpretation of energy flow,^{25–27} or the manner by which pressure across the scalae is focused at the organ of Corti complex (OCC) to supply energy to the traveling wave,^{28–31} whereas 1D models describe only the average pressure across the scalae.^{32,33}
The macromechanics of the cochlea are generally modeled as boundary value problems (BVPs), where the model equations involve partial differential equations (PDEs) without analytically known solutions, such as the NavierStokes equation. Such model equations can be tackled using numerics (e.g., the finite element method) or by making sufficient simplifying assumptions such that approximate analytic solutions can be found—the scala walls are rigid, the fluids are incompressible, etc.
Techniques based on the WentzelKramersBrillouin (WKB) approximation, also known as the LiouvilleGreen (LG) approximation,^{34–38} were introduced to the field of cochlear mechanics by Zweig and are among the most popular for achieving approximate, explicit closedform solutions for OCC motion, fluid pressure, and fluid velocity in a variety of cochlear models that match numerical solutions well across broad frequency and spatial ranges.^{4,5,9,21,28–31,39–41} Closedform explicit solutions allow for more easily interpreted model results. While exact solutions have been derived and studied for some cochlear models, e.g., implicit Green's function solutions^{17,18,42} or explicit Fourier transform solutions for box models,^{43} they are not as simple to qualitatively analyze.
As a nonexhaustive list, 1D and 2D WKB approximate solutions have offered: interpretations of limits on cochlear tuning,^{4} interpretations of intracochlear reflections and OAEs,^{10–13,15,31,44} interpretations of traveling wave mode coupling,^{22–24} and interpretations of the impact of active power generation in the cochlea.^{29,30} These approximations can also be applied to accelerate computations in more complex 3D cochlear models (e.g., Ref. 45) Robustness, ease of implementation, interpretability, and versatility have earned the WKB approximation its persistence in macromechanics modeling over half of a century.
With the passage of time, foundations of WKB techniques have largely disappeared from cochlear mechanics literature; as with any historical method, derivations, assumptions, implementation details, and the implications thereof have become implicit. This efficiency is useful for experienced readers but creates confusion for newer entrants to the field. In the case of WKB, not only are these objects often missing in contemporary literature, they are challenging to find in historic literature as well.
The relevance of the WKB approximation in cochlear models, historical and contemporary, owes it a foundational exposition. Fundamental understanding of the approximation can open the door to adaptations for probing particular questions with knowledge of its strengths and limitations. As such questions continue to arise with the publication of new data, especially with the advent of optical coherence tomography, this is all the more relevant.
Last, the recent passing of Egbert de Boer, Hendrikus Duifhuis, and Charles Steele—three pioneers of the application of the WKB approximation to cochlear mechanics—suggests a timeliness of such a presentation.
The essence of this tutorial is to present the fundamentals of WKB techniques in linear 1D and 2D cochlear mechanics models from an analytic perspective, covering derivations and details of implementation and performance.
I begin by describing general mathematical details of the WKB approximation agnostic to cochlear applications. This is followed by a description of the 1D and 2D BVPs for the box model of the cochlea. Derivations of the 1D and 2D WKB solutions to these BVPs follow. I then discuss the theory of the WKB traveling wave subspace (in terms of “WKB basis functions”) most often used in the study of intracochlear reflections and OAEs.
For readers interested in implementation rather than theory, Sec. VII discusses practical concerns. This includes discussion of several common methods for solving the dispersion relation for 2D box models, along with details of their performance across frequency and spatial ranges. This is followed by a comparison across methods and numerical solutions.
II. THE WKB APPROXIMATION
In this section, I will present the mathematical underpinnings of the WKB approximation. These abstract concepts will be applied to cochlear mechanics models specifically in Secs. III–VIII.
III. COCHLEAR MODEL
The WKB approximation can be applied to any cochlear model described by linear differential equations. In this tutorial, the focus is on one popular class of models—box and tapered box models—with geometry as displayed in Fig. 1. The model is derived by considering the cochlea as uncoiled and containing two scalae—scala vestibuli (SV) and scala tympani (ST).^{48,49} The scalae are separated by an infinitesimally thin plate, where the flexible OCC is the only portion of this plate capable of movement.
The cochlea's longitudinal axis (x) points toward the apex, terminating at the stapes at x = 0 and the helicotrema at x = L. The transverse axis (z) points toward SV with the OCC lying at z = 0. The crosssectional areas of SV and ST are equal to one another^{50} and vary along the longitudinal axis as S(x). The OCC width varies along the longitudinal axis as b(x). This model simplifies to the common box model when the scala walls are not curved and S(x) and b(x) are constant.
The model can be flattened to a 2D model as is represented geometrically in Fig. 1. This flattening amounts to representation of each quantity as its average over the radial dimension. It appears as a tapered box with height $ h ( x ) = S ( x ) / b ( x )$. Further flattening of the model to a 1D model involves representation of all quantities as being only dependent on x. This amounts to averaging quantities over transverse space.
A. Boundary conditions and assumptions
The boundary conditions are determined based on the following assumptions: (1) fluid does not flow in the normal direction toward or out of the scalae at the walls $ z = \xb1 h$ and helicotrema x = L, (2) the average pressure at x = 0 is some known constant P_{OW}, and (3) the OCC is mechanically described by a longitudinally varying point impedance, $ Z O C ( x )$ (or reciprocally as a point admittance $ Y O C = 1 / Z O C$).^{51} This quantity is complex and frequency dependent.
It should be noted that the modeled pressure and velocity will vary sinusoidally. Assuming linearity of the model,^{11,41,52–54} inputs at a given radian frequency, ω, will yield model responses at the same frequency. That is, all quantities will be of the form $ C ( x , z , \omega ) e j \omega t$. The timedependence is identical across all quantities and, therefore, it will generally be left implicit.
The fluid pressure is denoted by P(x,z), the longitudinal fluid velocity is denoted by $ u \u0307 ( x , z )$, and the transverse fluid velocity is denoted by $ w \u0307 ( x , z )$. The impedance describes the relationship between transmembrane pressure and transverse displacement at z = 0. Due to the symmetry of the model, transmembrane pressure $ p = 2 P$.
An assumption is also made regarding the character of the traveling wave solutions. With the input pressure appearing at the stapes, the traveling wave will primarily travel toward the apex. In this text, these will be referred to as apicaltraveling waves. At the helicotrema, reflection will occur and a wave traveling toward the stapes will be generated. Here, such waves will be referred to as basaltraveling waves.
For the most part, I will assume that basaltraveling waves are far smaller than apicaltraveling waves. In some models, this is achieved by letting $ L \u2192 \u221e$, in which case the WKB solutions will be identical to those arrived at in this tutorial. Basaltraveling waves will be considered in the context of WKB basis functions (Sec. VI), and this assumption will be removed.
B. Model equations
Another model equation can be derived for the average pressure in a cross section, i.e., for the 1D model in which pressure and velocity depend only on x. Over a small longitudinal cross section from x to $ x + \delta $, transverse fluid displacement occurs at a rate of approximately $ \delta b ( x ) w \u0307 ( x )$, as transverse fluid motion is generated only by the motion of the OCC, and $ \delta b ( x )$ is the approximate area of the OCC in this range.
IV. WKB SOLUTIONS FOR THE 1D MODEL
The Webster horn equation [Eq. (19)] was studied in the context of cochlear mechanics models as early as 1950.^{3} For general values of k(x), the PDE may not have a simple closedform analytic solution, but approximations such as constant k and simple forms for S can be used to yield explicit solutions.^{3} In the box model, where this simplifies to a wave/transmission line equation with varying wavenumber, an explicit solution is guaranteed but only in terms of retarded Green's functions.^{17,18,55,56} This motivates the development of a closedform, explicit approximate solution.
To arrive at approximate solutions, even further simplifying assumptions must be made. At first, these simplifications may appear as “hand waving” for the sake of mathematical convenience, but they will be kept track of and analyzed at the end of this section.
A. The zerothorder solution
B. The firstorder solution
Equations (34) and (38) are explicit formulas for pressure in terms of model parameters, geometry and wavenumber, found through density, frequency, and impedance [Eq. (20)]. Having solved for pressure, velocity of the OCC in the 1D model can be determined simply by dividing by the negative of the impedance.
C. The slowly varying parameter approximation
The term WKB assumption is often used to refer to the assumption that the wavenumber varies slowly in space relative to its own magnitude.^{4,5} However, the derivation above for first and secondorder WKB approximations never explicitly made this assumption. Where is the relationship between these two ideas?
With the dependence of k on the scala height and the impedance at the OCC [Eq. (20)], its rate of change is related to those of all of the model parameters. Thus,
WKB assumption: The parameters of the model vary slowly relative to their magnitudes.
When these conditions are not met, the WKB approximation breaks down. It is important to keep this assumption in mind when observing modeled responses.
Consideration of asymptotic behavior of the cochlea's traveling wave in light of Eq. (41) is instructive. At positions far basal to the best place, the response is said to be in the longwave region. Here, the wavenumber varies slowly in space, therefore, the lefthand term in Eq. (41) is very small.
At positions near the best place, the wavelength becomes smaller (wavenumber becomes larger) and changes more rapidly in space. This is known as the shortwave region. Here, the righthand term in Eq. (41) is very large. In balance, this assumption may be satisfied across a large portion of the frequency range.
On the other hand, reasonable smooth values for impedance will lead to quickly varying k around the region where stiffness and mass terms of the impedance cancel [see the dispersion relation of Eq. (20)]. In lossless cases, this leads to an infinite admittance and with small resistance still yields rapidly varying k.
D. The strong area assumption and the cochlear catastrophe
The presence of the wavenumber in this rewritten strong area assumption implies that there must be some balance between k and S to maintain the validity of the WKB approximation across the length of the cochlea. Reciprocally, appearance of the flattened scala height, h, in the formula for the wavenumber [Eq. (20)] implies that the derivative of scala area also impacts the WKB assumption for wavenumber variation in Eq. (41). Zweig and Shera discuss the implications of this balancing act between geometry and OCC impedance in detail and refer to the failure of models to account for this as the “cochlear catastrophe.”^{6} It is worth noting, however, that this catastrophe is only significant in the base in response to very lowfrequency stimuli, hence, box models without tapering reasonably satisfy the WKB approximation across most of space and frequency.
V. WKB SOLUTIONS FOR THE 2D MODEL
A. Pressure focusing
B. A higherorder 2D approximation
The above derivation arrives at Eq. (55) through two consecutive applications of the WKB approximation and neatly piggybacks off of 1D results for average pressure. However, this formula is not the only solution referred to in the literature as “the WKB solution” for a 2D model.
One derivation of this formula involves the solution for p as a formal power series approximation,^{21} inspired by the physics of surface waves.^{62} It is informative but lengthy, and an outline can be found in the Appendix. A second derivation of this formula follows from considering the EulerLagrange equations in a lossless box model (i.e., Z_{OC} purely imaginary).^{27} Neither derivation explicitly relies on the WKB approximation, although they do rely on the WKB assumption and feature the characteristic WKB phase term (the integral of the wavenumber). Intricate treatments of both derivations can be found online.^{78}
Although Eq. (56) behaves similarly to Eq. (55), its responses match numerical solutions better in the peak region. On the other hand, Eq. (56) only holds for box models where h is constant and does not allow for the modeling of cochlear tapering. Contemporary work is largely partial to the lowerorder approximation of Eq. 55.^{29–31,41,63} Differences in the behaviors between these solutions will be discussed in Sec. VIII.
C. Long and shortwave solutions
It is also instructive to consider the behavior of the solution in the longwave (small k, basal to best place/lower frequency than best frequency) and shortwave (large k, near best place/best frequency) limits.^{3,16,64,65} These approximations lie in the limiting behavior of the hyperbolic tangent for real argument $ a \u2208 \mathbb{R}$: $ tanh \u2009 a \u2248 a$ if a is small and $ tanh \u2009 a \u2248 1$ if a is large.
To visualize the differences between the longwave, shortwave, and WKB approximations, one can observe the effective height as k varies. Figure 2 shows h_{e} for the three solutions across various values of positive real k with h = 1 mm. It can be observed that the WKB solution for h_{e} exhibits a continuous switchoff between the long and shortwave approximations near the point where these solutions intersect. Behaviors of long and shortwave velocity responses are discussed in Sec. VIII.
VI. THE WKB TRAVELING WAVE SUBSPACE
Sections IV and V described derivations of explicit equations for pressure in 1D or 2D tapered box models via the WKB approximation. These formulas are valid where the model parameters do not change rapidly relative to their magnitudes and basaltraveling waves are negligible. However, WKB approximate solutions may not be easily derived for other cochlear models.
Numerical solutions, although more accurate to the dynamics, are generally challenging to interpret in comparison to WKB approximate solutions. This problem arises in, for example, the study of reflections in cochlear models—with only a numerical solution, how does one separate components of the response that are caused by incident waves from those resulting from reflected waves? This same breakdown may also be challenging in solving alternate cochlear models featuring, for example, nonlinearity. It is in this context that the theory of cochlear basis functions was developed.^{10}
Consider a solution to a 1D or 2D cochlear model that describes the pressure at the OCC in response to a singlefrequency stimulus. This is an infinitedimensional object, living in the Hilbert space, $H$, of smooth functions mapping from $ I = [ 0 , L ]$ (the interval of $\mathbb{R}$ along which the OCC is modeled to span) to $\u2102$.
Whereas the set of exact solutions is a subset of this infinitedimensional space, they will likely have qualitatively similar traveling wave forms for various perturbations to parameters and boundary conditions. Thus, they are likely to be wellapproximated as living in a lowerdimensional subspace containing functions that resemble cochlear traveling waves. This motivates the concept of a traveling wave subspace.
A. Theory of basis wave projection
In a pioneering work, Shera and Zweig introduce several sets of basis functions that generate a traveling wave subspace, including the WKB basis functions.^{10} I will develop the 1D box model WKB basis, but the method is just as well extended to other approximate solutions. This specification is for the sake of simplicity and because 1D box model WKB basis functions are the most commonly observed in the literature.^{10–15,44}
Of course, any exact solution to the BVP will not live in $W$, hence, the values for $ \psi \xb1$ found through this formula will not be constant. Thus, the projections are merely approximations that are best if the derivatives of $ \psi \xb1$ are sufficiently small.^{69}
Having these projection operators, one can formulate a numerical method for determining the apical and basaltraveling components of any pressure waveform by implementing derivatives as finite differences. The same process can be followed for other basis functions of approximate solutions, such as the shortwave solutions, longwave solutions, or WKB solutions in a tapered box model.
B. Applications to intracochlear reflections
Conversely, one can consider how basaltraveling waves reflect toward the apex via $ R \u2212$. With a passive cochlea driven at the stapes, this represents “reflections of reflections.” However, it is interesting to consider models where the cochlea is driven from a point along the length of the OCC ( $ x \u2260 0$).^{8,10,11} This could correspond to mechanical energy sources along the length of the OCC, present in the electromotile outer hair cells, which are likely responsible for many forms of OAEs.
Given that OAEs are measurable when the cochlea is driven at the stapes, there must be some significant portion of energy traveling toward the base without being entirely reflected. Some early modeling work on this topic predicted that the apical reflection coefficient is very large compared to the basal reflection coefficient ( $ R \u2212 \u226b R +$) such that basaltraveling energy would be significantly reflected before arriving back at the stapes.^{8,9} In this formulation, OAEs would have very low magnitudes. It was later argued by Shera and Zweig that the sizes of these quantities are highly dependent on the boundary conditions of the model^{10}—an important result to keep in mind for the modeling of OAEs.
C. Nonhomogeneous models and WKB solutions as a fundamental set
This formulation has broad applications in the modeling of intracochlear reflections and OAEs, where model equations can be manipulated into the form of Eq. (71) (see Refs. 11 and 12). In these cases, the forcing function, g, will generally represent sources of reflections such as random perturbations in impedance or nonlinearity. This interpretation is visible in Eq. (74), where g can be thought of as a kernel in the integral of the basis function traveling in the opposite direction of that for which it is a coefficient. That is, the size of the apicaltraveling component is modulated by the basaltraveling wave weighted by g and vice versa.
D. Example: Analytic treatment of roughness
There are various applications of WKB basis functions to the study of cochlear phenomena (several described in Refs. 11 and 12). In particular, values of g can be formulated to study sources of reflection, including nonlinear phenomena (e.g., distortion product OAEs). Here, I will provide a representative and important example—that of applying roughness to the cochlea's parameters.
VII. IMPLEMENTATION: SOLVING THE DISPERSION RELATION
In Secs. II–VI, I have described theoretical underpinnings for WKB solutions to 1D and 2D box and tapered box models. In this section, I discuss the challenges involved in implementation of the derived model equations in software.
At each position and frequency, solutions will exist for multiple values for k, but we will select only the most significant of such modes.^{72} As the velocity is loosely of the form $ e \u2212 jkx$, the solution should possess a positive real part to correspond to an apicaltraveling wave. As for the imaginary part, this leads to either dampening or amplification of the solution in x. Exponential growth is not physical in a passive cochlea, meaning that the imaginary part must be negative and the solution for k must lie in the fourth quadrant of the complex plane.
Moreover, of the solutions in this quadrant, the solution with the smallest (in magnitude) imaginary part is desired. A more negative imaginary part would lead to more severe exponential damping, therefore, the most significant solution has the least such damping.
In this section, I discuss the properties of the roots of f and the challenges that come in solving for physically reasonable roots. I, then, describe in detail three algorithms for finding k. The performance of these algorithms is discussed in Sec. VIII.
A. Root loci
Because the function f is continuous, a small variation of C should yield a small variation of the root position. Each continuous path traced by the roots with increasing x is called a root locus. With realistic impedance functions, the root loci form arcs in the fourth quadrant, traveling from the positive real axis to negative imaginary axis with increasing x.^{73} Fig. 3 shows four such root loci under one set of parameters, where each color corresponds to a different stimulus frequency and each circle is a root at a different x position (xresolution is $ 7 \u2009 \u2009 \mu m$). As x increases, the locus diverges from the real line and traverses clockwise toward the negative imaginary axis. At higher frequencies, the arc is broader and arrives at a larger negative imaginary value.
The WKB assumption is that this variation of k in x is slow such that tracing the continuous arc through the plane (possible with a fine enough resolution in x) would give the root of interest. However, with physically realistic parameters, one runs into multiple issues just past the peak region. In particular, near the location where stiffness and mass of Z_{OC} cancel [Eq. (82) and Table I], the admittance factor of C varies rapidly. Here, the WKB assumption breaks down, and a tracing of the root locus shows a rapid traversal of the arc near these positions. This can be observed in Fig. 3, where the roots appear sparse along the broad arc of the locus, indicating a much faster change in k than at the denser regions near the real and imaginary axes. In this region, insufficient resolution in x could not capture the continuous but rapid arc of the root locus and may, instead, yield convergence to a root in a different locus. This issue can be resolved either by uniformly refining resolution or refining resolution close to the resonant position.^{21}
Parameter .  Symbol .  Value . 

Mass  m  $ 1.5 \xd7 10 \u2212 3$ g/mm^{2} 
Resistance  r  $ 2 \xd7 10 \u2212 6$ Ns/mm^{3} 
Stiffness  s(x)  $ 10 e \u2212 0.2 x$ N/mm^{3} 
Scala height  h  1 mm 
Cochlea length  L  35 mm 
Fluid density  ρ  $ 10 \u2212 3$ g/mm^{3} 
Threshold on k finite difference  T  0.02 mm^{−2} 
Iterations of Newton's method (Algorithms 1 and 2)  M  20 
Iterations of contractive mappings (Algorithm 3)  M  20 
Points in xspace  N  1024 
Points in zspace (finite difference method)  N/A  16 
Parameter .  Symbol .  Value . 

Mass  m  $ 1.5 \xd7 10 \u2212 3$ g/mm^{2} 
Resistance  r  $ 2 \xd7 10 \u2212 6$ Ns/mm^{3} 
Stiffness  s(x)  $ 10 e \u2212 0.2 x$ N/mm^{3} 
Scala height  h  1 mm 
Cochlea length  L  35 mm 
Fluid density  ρ  $ 10 \u2212 3$ g/mm^{3} 
Threshold on k finite difference  T  0.02 mm^{−2} 
Iterations of Newton's method (Algorithms 1 and 2)  M  20 
Iterations of contractive mappings (Algorithm 3)  M  20 
Points in xspace  N  1024 
Points in zspace (finite difference method)  N/A  16 
B. Continuous longitudinally stepping algorithm
The goal is to begin by tracing a single root locus for f through the complex plane. Due to the number of possible roots at a given x, canonical rootfinding methods can cause trouble. Such methods require an intelligently chosen starting point so as not to converge to the wrong root or even a saddle point.
In this section, I will describe a class of algorithms for rootfinding that step across the longitudinal axis at each point, making an estimate for k informed by the estimate from the previous step.^{21,73} Here, x values are quantized into an Nlength vector with resolution $ \Delta x$. I will write the estimate for k at position x_{n} as $ k \u0302 n , \u2009 \u2009 n = 1 , 2 , \u2026 , N$. As the function f is itself x dependent, I will write $ f ( z ; x n )$ to refer to f at each position.
Starting at the very base, we are likely to be in the longwave region. This motivates the initial approximation of $ k \u0302 1 = k l w ( x 1 )$. This can be used as the initial value in a standard rootfinding algorithm such as NewtonRaphson or the Muller method, which are likely to converge to the correct root.
Stepping further along in x, the longwave approximation becomes poor. This indicates that we should not use this initial value forever. As in Fig. 3, wavenumbers within a single locus at subsequent x locations are likely close to one another—that is, $ k n \u2248 k n + 1$. The intuitive estimate is to use the solution at x_{n}, $ k \u0302 n$, as the starting point for the rootfinding method at $ x n + 1$ to find $ k \u0302 n + 1$.
This works so long as k is slowly varying, which is precisely the WKB assumption. However, there is a significant problem that occurs near the resonant point where stiffness and mass cancel, creating a rapid change in k (see the sparse regions of the arcs in Fig. 3).^{74} If we were to ignore this feature, we would simply follow the continuous root locus as in Fig. 3, tending toward solutions for k with large negative imaginary parts. This leads to falloff in the magnitude response that is far more rapid than what is observed in basilar membrane displacement data.^{27,75}
$ k \u0302 i c \u2190 k l w ( x 1 )$  $\u25b9$ Initialize using longwave k 
for $ n = 1 \u2192 N$ do $ z \u2190 h k \u0302 i c$  $\u25b9$ N is the number of steps in x space 
for $ m = 1 \u2192 M$ do $ z \u2190 z \u2212 f ( z ; x n ) f \u2032 ( z ; x n )$  $\u25b9$ M is the number of NewtonRaphson iterations 
end for  
$ k \u0302 n \u2190 z / h k \u0302 i c \u2190 k n$  $\u25b9$ Initial value for next step is current guess for k 
end for 
$ k \u0302 i c \u2190 k l w ( x 1 )$  $\u25b9$ Initialize using longwave k 
for $ n = 1 \u2192 N$ do $ z \u2190 h k \u0302 i c$  $\u25b9$ N is the number of steps in x space 
for $ m = 1 \u2192 M$ do $ z \u2190 z \u2212 f ( z ; x n ) f \u2032 ( z ; x n )$  $\u25b9$ M is the number of NewtonRaphson iterations 
end for  
$ k \u0302 n \u2190 z / h k \u0302 i c \u2190 k n$  $\u25b9$ Initial value for next step is current guess for k 
end for 
This rapid falloff past the peak is a result of a poor selection for k. In particular, the roots along the continuous locus do not correspond to dominant modes once their imaginary parts become sufficiently negative. Methods have been developed to counteract this problem by considering a continuous switchoff between dominance of two modes^{22–24} or discretely swapping the root locus being followed near the resonant position.^{20,21} An example of the latter type is discussed in Sec. VII C.
C. Discontinuous longitudinally stepping algorithm
To account for the issues found in tracing the continuous root locus near the resonant position, one can introduce a discontinuity into the xstepping algorithm. This is performed by changing the initial value for the rootfinding algorithm near the resonant position, where the WKB assumption breaks down. To do so, we seek a better guess for a root near this position.
The discontinuous xstepping algorithm traces the continuous root locus up to some x_{d} at which it is determined that the WKB assumption is being violated. This can be determined before simulation^{21} or on the fly by observing the rate of change of the wavenumber (in discrete space, the finite difference) at each step. When the WKB assumption holds, this rate should be small. Picking some threshold T > 0, the xstepping method is paused once $  k \u0302 n \u2212 k \u0302 n \u2212 1  / \Delta x > T$.
After this point, k_{d} of Eq. (81) is used as an initial step in the rootfinding algorithm. If $  k \u0302 n \u2212 k \u0302 n \u2212 1  / \Delta x > T$ still holds, then the WKB assumption is violated, and the pressure at this position is set equal to its value at the last computed position ( $ p n = p n \u2212 1$). At each subsequent step, it is determined whether $ k \u0302$ satisfies this threshold—it will eventually do so, at which point we continue the locustracing process along this second locus. Pseudocode for this algorithm is presented in Algorithm 2.
$ k \u0302 i c \u2190 k l w ( x 1 )$  $\u25b9$ Initialize using longwave k 
for $ n = 1 \u2192 N$ do  $\u25b9$ N is the number of steps in x space 
$ z \u2190 h k \u0302 i c$  
for $ m = 1 \u2192 M$ do $ z \u2190 z \u2212 f ( z ) f \u2032 ( z )$  $\u25b9$ M is the number of NewtonRaphson iterations 
end for  
if $  z / h \u2212 k i c  \Delta x < T$ then  $\u25b9$ Check for validity of WKB condition 
$ k \u0302 n \u2190 z / h$  $\u25b9$ Initial value for next step is current guess for k 
$ k \u0302 i c \u2190 k n$  
else  
$ k \u0302 n \u2190 NaN$  $\u25b9$ Pressure and velocity should not be computed here 
$ k \u0302 i c \u2190 k d$  $\u25b9$ Guess for k after the discontinuity 
end if  
end for 
$ k \u0302 i c \u2190 k l w ( x 1 )$  $\u25b9$ Initialize using longwave k 
for $ n = 1 \u2192 N$ do  $\u25b9$ N is the number of steps in x space 
$ z \u2190 h k \u0302 i c$  
for $ m = 1 \u2192 M$ do $ z \u2190 z \u2212 f ( z ) f \u2032 ( z )$  $\u25b9$ M is the number of NewtonRaphson iterations 
end for  
if $  z / h \u2212 k i c  \Delta x < T$ then  $\u25b9$ Check for validity of WKB condition 
$ k \u0302 n \u2190 z / h$  $\u25b9$ Initial value for next step is current guess for k 
$ k \u0302 i c \u2190 k n$  
else  
$ k \u0302 n \u2190 NaN$  $\u25b9$ Pressure and velocity should not be computed here 
$ k \u0302 i c \u2190 k d$  $\u25b9$ Guess for k after the discontinuity 
end if  
end for 
D. The fixed point algorithm
One alternative to the longitudinally stepping class of algorithms is a fixed point algorithm in which two distinct relationships between k and α (the pressure focusing factor) are used.^{13,63} The first such relationship is Eq. (52), which gives α in terms of k. The second relationship, giving k in terms of α, is Eq. (23). Any valid k value must satisfy both equations.
Fixed point methods are based on the Banach fixed point theorem,^{77} which states that repeated application of a contractive function, g, will converge to a fixed point of said function, i.e., a point where g(x) = x. Mathematical details are omitted here for the sake of brevity.
The fixed point method for this problem works by starting with the longwave approximation at every frequencylocation pair, $ k \u0302 = k l w$. Then, $ k \u0302$ is plugged into Eq. (52) to find an the approximate pressure focusing factor, $ \alpha \u0302$, and then $ \alpha \u0302$ is plugged into Eq. (23) to find a new wavenumber approximation, $ k \u0302$. This is repeated for some number of iterations. Pseudocode for this algorithm is shown in Algorithm 3.
$ k \u0302 \u2190 k l w$  $\u25b9$ Here, $ k \u0302$ is a vector with an index for each position 
for $ m = 1 \u2192 M$ do  $\u25b9$ M is the number of fixed point iterations 
$ \alpha \u0302 \u2190 h k \u0302 tanh \u2009 k h$  $\u25b9$ Pressure focusing vector update [Eq. (52)] 
$ k \u0302 \u2190 \u2212 2 j \omega \rho \alpha \u0302 h Z O C$  $\u25b9$ Wavenumber update [Eq. (23)] 
if $ R [ k \u0302 ] < 0$ then  $\u25b9$ Ensure the root is for an apicaltraveling wave ( $R$ gives the real part) 
$ k \u0302 \u2190 \u2212 k \u0302$  
end if  
end for 
$ k \u0302 \u2190 k l w$  $\u25b9$ Here, $ k \u0302$ is a vector with an index for each position 
for $ m = 1 \u2192 M$ do  $\u25b9$ M is the number of fixed point iterations 
$ \alpha \u0302 \u2190 h k \u0302 tanh \u2009 k h$  $\u25b9$ Pressure focusing vector update [Eq. (52)] 
$ k \u0302 \u2190 \u2212 2 j \omega \rho \alpha \u0302 h Z O C$  $\u25b9$ Wavenumber update [Eq. (23)] 
if $ R [ k \u0302 ] < 0$ then  $\u25b9$ Ensure the root is for an apicaltraveling wave ( $R$ gives the real part) 
$ k \u0302 \u2190 \u2212 k \u0302$  
end if  
end for 
This method works under the assumption that it converges to the correct value of k, which depends on the properties of the mappings between α and k. If the mappings are not (at least locally) contractive, then no convergence is guaranteed. On the other hand, if there are multiple fixed points, certain choices of initial conditions may lead to convergence to an undesired k. Performance of these three algorithms will be discussed in Sec. VIII.
VIII. BEHAVIOR OF WKB APPROXIMATE SOLUTIONS
A. WKB solutions for the 1D box model
In Sec. IV, I derived WKB solutions to the 1D BVP up to the zeroth [Eq. (34)] and first [Eq. (38)] orders. In Fig. 4, these solutions are shown in response to a 2 kHz stimulus frequency and compared to numerical results.
It can be observed that the zerothorder approximation overestimates the magnitude of the response near the peak and exhibits more phase accumulation than the numerical solution. On the other hand, the firstorder WKB approximation matches the numerical solution remarkably well across space in phase and magnitude. The two orders of solution differ only by a factor of $ k$, which is real valued for small x, explaining the similarity in phase.
B. Longwave and shortwave solutions
To contextualize findings for the 2D WKB solutions, it is useful to observe the performance of the long and shortwave approximate solutions to the 2D box model (see Sec. V C as well as Refs. 3, 16, and 64). These solutions are valid for regions of small real k and large real k, respectively, but as k is complex valued and varies nonmonotonically across space/frequency (see Fig. 3), it is not immediately clear in which regions these approximations will best match numerical solutions.
Figure 5 shows the longwave and shortwave solutions to the 2D box model alongside a numerical solution. The longwave response matches the numerical solution well at more basal positions, where the wavenumber is small and real (Fig. 3), but poorly matches the numerical solutions near or above the peak.
The shortwave solution matches the numerical solution well only in a small spatial range near the peak. Neither approximation matches the numerical solution past the peak where the magnitude falloff in the numerical solution becomes slower. This region is termed the cutoff region.^{22} In the context of the root loci (Fig. 3), neither approximation should be expected to hold well in the cutoff region, where the roots approach the negative imaginary axis, as the asymptotic forms of the hyperbolic tangent used in their derivations are only valid for real argument. The shortwave solution exhibiting a sign change in its group velocity just past the peak, is one dramatic consequence of this breakdown.
C. Performance of wavenumberfinding algorithms
Before comparing 2D WKB approximations to numerical results, it is first important to assess the methods for determining the wavenumber k in the 2D case. This is performed by observing velocity responses at the OCC derived from the 2D WKB approximation of Eq. (56), using three methods for finding the wavenumber: (1) Algorithm 1, an xstepping algorithm that does not account for the discontinuity, amounting to following a root locus as in Fig. 3; (2) Algorithm 2, an xstepping algorithm that does account for discontinuity, via thresholding the finite difference as described above; and (3) Algorithm 3, the fixed point method.
Figure 6 contrasts the xstepping methods depending on whether discontinuity is accounted for. The velocity responses show identical behavior up to a position slightly past the peak, where the finite difference in k becomes sufficiently large such that a discontinuity is registered. After this point, the falloff in velocity amplitude is slower than that if the discontinuity were not considered. This slower falloff is observed in the cutoff region of numerical results, suggesting that the discontinuous method yields a more reasonable choice for k past the peak.^{20,21,27,73} Comparison to numerics is presented in Sec. VIII D.
Recall that these algorithms are designed to solve a rootfinding problem for function f of Eq. (76). The root loci for f in the upperright panel of Fig. 6 show that for the discontinuous method, the traversal of the locus is halted as the root pattern begins to appear sparser (i.e., faster change in k). As described in Sec. VII, the discontinuous algorithm then assumes a small negative imaginary root (transition shown by the dashed blue arrow), which yields less rapid falloff in the cutoff region than the larger negative imaginary component found by following the locus continuously.
The bottomright panel of Fig. 6 serves to show that the two methods are correctly converging to roots of f at each given x. Using both methods, the value of f(kh) is less than $ 10 \u2212 10$ in magnitude at all x—this stresses the fact that not all roots lie on the same continuous locus.
Figure 7 shows these same results but now alongside the results obtained via the fixed point method (Algorithm 3). These results show similar behavior in velocity magnitude to the continuous xstepping solution, but the phase accumulates more cycles.
Observation of the root locus and f(kh) for the fixed point algorithm reveals strange behavior in the cutoff region. While the fixed point method's root locus follows that of the xstepping method for a large range of x, it erratically jumps around the complex plane (including to the third quadrant) past the peak. This corresponds to a nonzero value of f(kh) at these positions as well (see the bottomright panel of Fig. 7), showing that Algorithm 3 has not correctly converged to a root of the function.
This is anecdotal justification of the validity of this method in the longwave and shortwave regions but not in the cutoff region—a drawback of the fixed point method. The failure of the fixed point algorithm to converge in the cutoff region has been noted before, e.g., in Appendix D of Ref. 63. Still, because of the relative speed of this method's convergence and the fact that it does not need a fine resolution for x, it has experienced use in many modern works where performance in the cutoff reason is not critical to model results.^{31,41,63}
D. WKB solutions for the 2D box model
In Sec. VIII C, it was shown that the wavenumberfinding algorithm that accounts for discontinuities in k converges to roots across space and yields velocity responses that qualitatively resemble numerical results past the peak. This informs the choice to use this algorithm for comparison to numerical results.
In Sec. V, two WKB approximate solutions were presented—the lowerorder solution of Eq. (55) and the higherorder solution of Eq. (56). Figure 8 shows solutions according to both of these equations alongside numerical solutions to the 2D BVP.
Both approximate solutions resemble the numerical solutions across space, including at positions past the peak where the slope of the response rapidly changes. This contrasts with the approximate solutions derived from the two alternate rootfinding methods, as observed in Fig. 7.
The lowerorder solution slightly overestimates tuning at the peak as compared to the higherorder solution, but the solutions are otherwise nearly identical. They differ most significantly from numerical solutions in their phase responses. Although they are characteristically similar, both approximate solutions lead the numerical solutions by about 0.1 cycles at apical positions where the phase varies slowly.
It is important to note that these results are sensitive to the choice of threshold, T, in Algorithm 2. A large threshold leads to a registration of a discontinuity at a more apical position. This means that the slopes of the magnitude and phase will change at a more apical location than in the numerical solution. Similarly, a lower value of T will move the discontinuity further basal. For reasonable values of T, WKB solutions will be qualitatively similar to numerical solutions in magnitude and phase slopes, but they may differ quantitatively as a result of the shift in the position at which the discontinuity is registered.
IX. CONCLUSIONS
The WKB approximation provides compact closedform approximate solutions for 1D and 2D cochlear macromechanic models that match numerical solutions well within the entire region of response (all x's and ω's) except a small region near the resonant position/frequency. These solutions are easily implemented and interpreted and allow qualitative and quantitative insight into the manner by which physical parameters (impedance variations, scala area, etc.) alter the apicaltraveling wave. Through the method of WKB wavespace projection, the solutions also facilitate further interpretation of modeled responses as a superposition of apical and basaltraveling waves, allowing for the study of intracochlear reflections and OAEs.
The forms of the WKB approximate solutions for cochlear box models were developed decades ago, thus, one may reasonably ask: what is the importance of WKB approximate solutions in contemporary times? In fact, many important insights have been gleamed from WKB solutions in recent literature. Below are just a few such contributions from the past three years.
The inclusion of spatially varying scala dimensions to a 2D box model as that in Eq. (55) has been shown by Altoè and Shera to be important for the achievement of substantial OCC velocity at the apex in response to stimuli at the base.^{63} They analyzed how tapering scala height could introduce an amplification factor that boosts responses at the apex relative to a uniformheight model, resolving losses that occur in the traveling wave as it makes its way to the apex.
Recent micromechanical findings made through optical coherence tomography have also inspired applications of the WKB solutions to the ostensibly macromechanical box model. In particular, motion at the outer hair cellDeiters cell junction in the organ of Corti appears to move about 90° out of phase with basilar membrane motion within the same longitudinal cross section. Implementing this as a modification to the impedance term, Altoè and Shera have used the WKB solutions as derived in this tutorial to model the impact of such a phenomenon.^{41} They arrive at an alternative interpretation of cochlear amplification, in which power may be supplied to the fluid rather than directly to the basilar membrane.
Recent work by Sisto et al. used the WKB approximation in studying the level dependence of the OCC admittance, assumed to arise from outer hair cell motility.^{29,30} Paying special attention to (a) the pressure focusing phenomenon described above and (b) the viscosity at the OCC–fluid interface, they have found that substantially leveldependent admittance is not required to obtain the impressive dynamic range of the cochlea.
With much still to learn about the mechanics of the cochlea, the powerful analytic tool offered by the WKB approximation is among the strongest we have because of its interpretability, versatility, and simplicity of computation. With the foundations discussed in this tutorial, derivations and implementation details can be modified to tackle contemporary questions as they continue to arise.
ACKNOWLEDGMENTS
I would like to thank Dr. Elizabeth S. Olson and the JASA reviewers for providing edits and comments for this tutorial, as well as National Institute on Deafness and Other Communication Disorders (F31 DC02062102).
AUTHOR DECLARATIONS
Conflict of Interest
The author has no conflicts to disclose.
Data Availability
No data were used in this tutorial article. Scripts used in generating the presented figures are available from the author upon reasonable request.
APPENDIX: THE HIGHERORDER 2D MODEL—A SERIES SOLUTION APPROACH
In this appendix, I provide a derivation of the higherorder “WKB solution” of Eq. (56). The following approach is modified from that of Viergever in his 1980 book, Mechanics of the Inner Ear: A Mathematical Approach.^{21} It relies on a transformation of the coordinates of the pressure BVP and subsequent application of a WKBadjacent ansatz (but not precisely the WKB method).
The method consists of the following steps:

Change the variables of the BVP in pressure such that terms relating to the model parameters appear in the PDE rather than only in the boundary conditions; and
 write a form for the solution to this new PDE as$ A ( x , \zeta ) cosh [ \kappa ( x ) ( H \u2212 \zeta ) ] e jKg ( x ) .$
That is, assume that the z (here, reparameterized as ζ) contribution is hyperbolic and there is a wave in x. The product with arbitrary A(x,y) means this is performed without loss of generality.

Assume a series solution for A and plug it into the ODE to obtain a system of PDEs; and

solve for A up to first order, plug back into the ansatz, and undo the change of variables to solve for pressure.
A detailed outline is presented below, but certain steps feature highly nontrivial computations. Full exposition of these computations is available online.^{78}
1. Setting up the BVP
The method followed in this section relies on multiple changes of variables and definitions of new parameters. As such, it can be difficult to keep straight the meanings and units of the various variables and parameters at play. Table II serves as a reference for the objects introduced in the derivation.
Symbol .  Significance .  Units . 

Z_{0}  Arbitrary reference impedance used to the simplify series solution  Pa s/mm 
K  Reference wavenumber used to simplify the series solution  1/mm 
$ f 2 ( x )$  $ Z 0 / Z O C ( x )$, used so simplify xdependence of the PDE  Unitless 
ζ  Kz, Nondimensionalized transverse coordinate  Unitless 
H  Kh, Nondimensionalized scala height  Unitless 
$ Q ( x , \zeta )$  Pressure written in terms of the nondimensionalized transverse coordinate  Pa 
$ A ( x , \zeta )$  Auxiliary pressure variable that controls the magnitude of pressure at the OCC, to be solved for in the simplified BVP  Pa 
$ \kappa ( x )$  Controls the xdependence of transverse pressure variations, to be solved for in the simplified BVP  Unitless 
g(x)  Controls the wavenumber of the traveling wave, to be solved for in the simplified BVP  mm 
Symbol .  Significance .  Units . 

Z_{0}  Arbitrary reference impedance used to the simplify series solution  Pa s/mm 
K  Reference wavenumber used to simplify the series solution  1/mm 
$ f 2 ( x )$  $ Z 0 / Z O C ( x )$, used so simplify xdependence of the PDE  Unitless 
ζ  Kz, Nondimensionalized transverse coordinate  Unitless 
H  Kh, Nondimensionalized scala height  Unitless 
$ Q ( x , \zeta )$  Pressure written in terms of the nondimensionalized transverse coordinate  Pa 
$ A ( x , \zeta )$  Auxiliary pressure variable that controls the magnitude of pressure at the OCC, to be solved for in the simplified BVP  Pa 
$ \kappa ( x )$  Controls the xdependence of transverse pressure variations, to be solved for in the simplified BVP  Unitless 
g(x)  Controls the wavenumber of the traveling wave, to be solved for in the simplified BVP  mm 
The exponential suggests a traveling wave in x, where A modulates the amplitude of this wave. However, this is not actually an assumption of a wave solution—as A, g, and κ are unknown functions of x and ζ, any function can be represented in this fashion without loss of generality.
One might wonder why we have chosen to introduce so many new terms into these equations. Although this may initially seem to complicate the BVP, it eventually leads to the most mathematically tractable solution method.
Alongside Table II, it may help to “look into the future” to see what these newly defined variables will become. The variable κ will be found to be the nondimensionalized wavenumber and g will be found to be the integral of the wavenumber. The free parameter K will eventually be the variable of our formal power series [similar to δ from Eq. (2)].
Solving this ODE in the auxiliary pressure, A, is the new goal. With a solution for A, we can find Q and, finally, P.
2. A series solution for auxiliary pressure
3. Finding a first approximation for pressure
Theoretically, this facilitates solution for A_{n} for any n as well. On the other hand, the series approximation concludes that the higher n terms should be small if K is large relative to its own rate of change (analogous to the WKB assumption).
The reference constant, K, was defined as $ K 2 = \u2212 2 j \rho \omega / h Z 0$, where Z_{0} was a second reference constant such that $ f 2 = Z 0 Y O C$. Because Z_{0} was entirely arbitrary, I am free to choose $ Z 0 = \u2212 2 j \rho \omega h \u2212 1$ such that K = 1 mm^{−1}. This also gives H = 1 mm^{−1} × h [unitless], ζ = 1 mm^{−1} × z [unitless] and Q = P (Pa).