The Wentzel-Kramers-Brillouin (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 closed-form 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 two-dimensional cochlear macromechanics. This includes mathematical background, rigorous derivation and details of its implementation in software.

Models of one- (1-D) and two-dimensional (2-D) macromechanics have offered some of the most significant interpretations of cochlear physiology, historic and modern. It is intuitive that three-dimensional (3-D) models should offer more physically realistic results than 2-D or 1-D 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 1-D models include the existence and character of stapes-driven 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 long-wave region to in vivo cochlear responses make 1-D 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 1-D models than in 2-D or 3-D models.10–13 

On the other hand, 2-D models allow for more physical results in the short-wave and cutoff regions of the cochlear response than 1-D models,16–21 allowing for more complete exploration of potential mode-coupling phenomena in the traveling wave.22–24 Moreover, 2-D 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 1-D 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 Navier-Stokes 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 Wentzel-Kramers-Brillouin (WKB) approximation, also known as the Liouville-Green (LG) approximation,34–38 were introduced to the field of cochlear mechanics by Zweig and are among the most popular for achieving approximate, explicit closed-form 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 Closed-form 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 solutions17,18,42 or explicit Fourier transform solutions for box models,43 they are not as simple to qualitatively analyze.

As a nonexhaustive list, 1-D and 2-D 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 3-D 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 1-D and 2-D 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 1-D and 2-D BVPs for the box model of the cochlea. Derivations of the 1-D and 2-D 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 2-D box models, along with details of their performance across frequency and spatial ranges. This is followed by a comparison across methods and numerical solutions.

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.

Consider a homogeneous linear nth-order ordinary differential equation (ODE) of the form
ϵ d n y d x + a n 1 ( x ) d n 1 y d x n 1 + + a 1 ( x ) d y d x + a 0 ( x ) y = 0 ,
(1)
where the coefficient functions ai, i = 1 , 2 , , n 1, are n-times continuously differentiable functions on some interval I , and ϵ is presumed to be small relative to the magnitudes of the other coefficient functions. The coefficient functions may be complex valued.
Consider an ansatz for a solution to Eq. (1) as the exponential of a formal power series in δ ,
y ( x ) = exp [ 1 δ m = 0 δ m C m ( x ) ] ,
(2)
where Cm, m = 0 , 1 , 2 , , are n-times continuously differentiable functions on I, and it is assumed that the series can be differentiated term-wise.34,35,37,38,46 For the series to converge,47  δ should be small and Cm and their derivatives must fall off exponentially in magnitude across the real line with increasing m. That is,
| δ m C m + 1 ( x ) | | δ m 1 C m ( x ) | , m = 0 , 1 , 2 , .
(3)
The ansatz, when plugged into Eq. (1), yields a differential equation containing infinitely many unknowns, Cm. The Mth-order WKB approximation is made by truncating this series up to the Mth term. This is valid so long as all terms at indices higher than M are much smaller than one on I. That is,
| δ M C M + 1 ( x ) | 1 .
(4)
It is common that the first-order WKB approximation is simply called “the WKB approximate solution.”36 

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.

FIG. 1.

(Color online) (A) Geometry of the 3-D box model, having length L, scala height h, and width b is shown. The longitudinal, radial, and transverse directions are x, y, and z, respectively. (B) Geometry of the corresponding 2-D box model is displayed. In the tapered box model, h may be a function of x. SV, scala vestibuli; ST, scala tympani; OW, oval window.

FIG. 1.

(Color online) (A) Geometry of the 3-D box model, having length L, scala height h, and width b is shown. The longitudinal, radial, and transverse directions are x, y, and z, respectively. (B) Geometry of the corresponding 2-D box model is displayed. In the tapered box model, h may be a function of x. SV, scala vestibuli; ST, scala tympani; OW, oval window.

Close modal

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 cross-sectional areas of SV and ST are equal to one another50 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 2-D 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 1-D model involves representation of all quantities as being only dependent on x. This amounts to averaging quantities over transverse space.

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 = ± h and helicotrema x = L, (2) the average pressure at x = 0 is some known constant POW, 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 , ω ) e j ω t. The time-dependence 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 ̇ ( x , z ), and the transverse fluid velocity is denoted by w ̇ ( 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.

In the 2-D tapered box model, the boundary conditions can be written as
1 h ( 0 ) 0 h ( 0 ) P ( 0 , z ) d z = P O W ,
(5)
P x ( L , z ) = 0 ,
(6)
P z ( x , h ( x ) ) = 0 ,
(7)
P ( x , 0 ) = Z O C ( x ) 2 w ̇ ( x , 0 ) ,
(8)
where the negative sign in the impedance relation comes from positive pressure in SV applying a –z-direction (downward) force on the OCC. These boundary conditions are equally valid in the 1-D model simply by ignoring dependence on z.

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 apical-traveling 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 basal-traveling waves.

For the most part, I will assume that basal-traveling waves are far smaller than apical-traveling waves. In some models, this is achieved by letting L , in which case the WKB solutions will be identical to those arrived at in this tutorial. Basal-traveling waves will be considered in the context of WKB basis functions (Sec. VI), and this assumption will be removed.

With the geometry described, the model equations can now be developed. The fluid in the scalae is modeled as incompressible, irrotational, linear, and inviscid. The fluid velocity vector, v = ( u ̇ w ̇ ) T, in a volume satisfies the continuity equation, which is given by
ρ t + ( ρ v ) = 0 ,
(9)
where ρ is the fluid density. This represents that within a differential volume, the change in fluid mass in the region is accompanied by an equal and opposite divergence of that fluid into/out of the region.
In an incompressible fluid, the mass of the fluid (and, thereby ρ) in any region is constant, simplifying the equation to
v = 0.
(10)
An irrotational field is also a conservative field, thus, the velocity field can be written as the gradient of some scalar field, ϕ. This velocity potential is thereby defined by
ϕ = v .
(11)
Taking the divergence of both sides and applying Eq. (10) yields the Laplace equation in velocity potential:
2 ϕ = 0.
(12)
The Navier-Stokes equation in an inviscid, incompressible, linear, and irrotational fluid is
ρ v t + P = ρ ϕ t + P = 0.
(13)
This gives
P = ρ ϕ ̇ ,
(14)
where the overhead dot indicates a partial derivative in time. Taking the Laplacian of both sides and recalling that ϕ satisfies the Laplace equation, we arrive at a Laplace equation in P such that
2 P = 0.
(15)
The Laplace equation then also holds for transmembrane pressure, p = 2 P.

Another model equation can be derived for the average pressure in a cross section, i.e., for the 1-D model in which pressure and velocity depend only on x. Over a small longitudinal cross section from x to x + δ, transverse fluid displacement occurs at a rate of approximately δ b ( x ) w ̇ ( x ), as transverse fluid motion is generated only by the motion of the OCC, and δ b ( x ) is the approximate area of the OCC in this range.

Longitudinally, fluid enters the region at rate S ( x ) u ̇ ( x ) and exits at rate S ( x + δ ) u ̇ ( x + δ ). Exiting fluid must be the sum of entering fluid and fluid displaced by transverse OCC motion:
S ( x + δ ) u ̇ ( x + δ ) = S ( x ) u ̇ ( x ) + δ b ( x ) w ̇ ( x ) .
This can be manipulated into the form of a difference quotient,
S ( x + δ ) u ̇ ( x + δ ) S ( x ) u ̇ ( x ) δ = b ( x ) w ̇ ( x ) .
Letting δ 0, the left-hand side is recognized as a derivative in x:
x [ S u ̇ ] = b w ̇ .
(16)
This identity can be used to arrive at a differential equation in pressure. This begins with simplifying the Navier-Stokes equation [Eq. (13)] to a 1-D equation, multiplying it by the cross-sectional area and differentiating in x such that
x [ S P x ] + ρ x [ S u ̇ t ] = 0 ,
(17)
where I have used the fact that u ̇ = ϕ / x.
Replacing the time derivative by-product with j ω and applying Eq. (16) results in
x [ S P x ] + j ω ρ b w ̇ = 0.
(18)
To write this entirely in terms of transmembrane pressure, p, I can use p = 2 P and the 1-D model's point-impedance boundary condition, p = Z O C w ̇. Doing so and dividing by S gives the Webster horn equation for the 1-D model:
1 S x [ S p x ] + k 2 p = 0 ,
(19)
k 2 ( x ) = 2 j ω ρ Z O C ( x ) h ( x ) .
(20)
This derivation can be readily modified to apply to the 2-D model as well—one replaces the 1-D model's longitudinal velocity, u ̇, and pressure, P, with those of the 2-D model averaged across the transverse dimension and replaces the 1-D model's transverse velocity, w ̇, with that of the 2-D model at z = 0. The derivation holds identically until the final step. In the 2-D model, the boundary condition is that transverse OCC velocity is related to the pressure at z = 0, not the average pressure. Writing the average pressure as P ¯ (or p ¯ = 2 P ¯), the pressure focusing factor is defined as
α ( x ) = p ( x , 0 ) p ¯ ( x ) ,
(21)
which is the ratio between the pressure focused at the OCC and average transmembrane pressure in the cross section. The 2-D model Webster horn equation is then
1 S x [ S p ¯ x ] + k 2 D 2 p ¯ = 0 ,
(22)
k 2 D 2 ( x ) = 2 j ω ρ α ( x ) Z O C ( x ) h ( x ) .
(23)

In the box model where S is constant, Eqs. (19) and 22 degenerate to wave equations with variable wavenumbers, and Eqs. (20) and 23 are dispersion relations.

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 closed-form 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 closed-form, explicit approximate solution.

The WKB approximate solution for the horn equation57–59 can be considered by placing the equation into the standard form of a linear differential equation:
2 p x 2 + S S p x + k 2 p = 0 ,
(24)
where “.′” denotes the spatial derivative. Comparing with Eq. (1), we have ϵ = 1.
Putting in the WKB ansatz [Eq. (2)] with δ = ϵ = 1 gives
exp ( m = 0 C m ) [ m = 0 C m + ( m = 0 C m ) 2 + S S m = 0 C m + k 2 ] = 0.
(25)
As the exponential term is never zero and leads every term, I can divide through by it. I choose to keep only terms involving C0 and C1. By the asymptotic assumptions of the WKB approximation, C 1 C 0, the terms should decrease as further derivatives are taken such that C 1 and ( C 1 ) 2 are negligible in comparison to lower-order terms. This results in
C 0 + ( C 0 ) 2 + 2 C 0 C 1 + S S ( C 0 + C 1 ) = k 2 .
(26)

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.

The most stringent approximation made in these derivations is the strong area assumption,
| S S C 0 | | C 0 | 2 ,
(27)
an assumption about the physical parameters of the system that can be interrogated once a first approximation of C0 has been found. This is trivially true in a box model no matter what the value of C0 is, as in this case S = 0.
A consequence of this strong area assumption is the weak area assumption,
| S S C 1 | | S S C 0 | | C 0 | 2 ,
(28)
which is justified by the asymptotic assumptions on the WKB series. Application of this weak area assumption to Eq. (26) gives the first-order WKB ODE:
C 0 + ( C 0 ) 2 + 2 C 0 C 1 + S S C 0 = k 2 .
(29)
To obtain the zeroth-order WKB ODE, I perform a further reduction based on the same asymptotic decay assumptions. First, I make the second derivative assumption,
| C 0 | | C 0 | 2 ,
(30)
such that I am justified in ignoring the first summand of Eq. (29). Also ignoring the remaining first-order term of 2 C 0 C 1 and applying the strong area assumption, Eq. (29) reduces to
( C 0 ) 2 = k 2 .
(31)
With Eqs. (31) and 29, zeroth- and first-order WKB approximations for the 1-D model can be found.
The ODE in Eq. (31) is quickly solved by taking the square root of both sides and integrating:
C 0 = ± j 0 x k ( ξ ) d ξ .
(32)
Plugging into the WKB ansatz, the zeroth-order solution is
p 0 ( x ) = A e j 0 x k ( ξ ) d ξ + B e j 0 x k ( ξ ) d ξ , A , B .
(33)
This is recognized as a superposition of two traveling waves with the first summand traveling toward the apex and the second summand traveling toward the base. As stated in Sec. III A, the basal-traveling waves caused by reflection at the helicotrema are modeled as being negligible, thus, B = 0.
Application of the boundary condition at the oval window [Eq. (5)] with B = 0 gives A = P O W. The zeroth-order solution is then
p 0 ( x ) = P O W e j 0 x k ( ξ ) d ξ .
(34)
Putting in the value for C0 found in Eq. (32) to the first-order WKB approximate ODE of Eq. (29), I have
± j k k 2  ±  2 j k C 1  ±  S S j k = k 2 .
Solving for C 1 results in
C 1 = k 2 k S 2 S ,
(35)
and integrating both sides,
C 1 = 1 2 log ( S k ) .
(36)
Plugging the found values of C0 and C1 into the WKB ansatz gives the first-order WKB approximate solution,
p 1 ( x ) = A S ( x ) k ( x ) e j 0 x k ( ξ ) d ξ + B S ( x ) k ( x ) e j 0 x k ( ξ ) d ξ , A , B .
(37)
Once again, B = 0 by the assumption that basal-traveling waves are negligible. A is found by applying the boundary condition at the oval window and writing S ( 0 ) = S 0 , k ( 0 ) = k 0 such that
p 1 ( x ) = P O W S 0 k 0 S ( x ) k ( x ) e j 0 x k ( ξ ) d ξ .
(38)
In the box model, as S ( 0 ) = S ( x ) for all x, the ratio inside the square root simplifies to k 0 / k ( x ).

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 1-D model can be determined simply by dividing by the negative of the impedance.

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 second-order WKB approximations never explicitly made this assumption. Where is the relationship between these two ideas?

For the WKB method to be valid, the terms Cn in the series must decrease monotonically. In particular, | C 1 | | C 0 |. Using the box model case for the sake of simplicity (no dependence on S), consider this relationship with the values of C0 and C1 from the 1-D model derived above,
| 1 2 ln k | | j 0 x k ( ξ ) d ξ | .
(39)
The left-hand expression can also be written as an integral from zero to x, and pulling out the constant-modulus factors gives
1 2 | 0 x k ( ξ ) k ( ξ ) d ξ | | 0 x k ( ξ ) d ξ | .
(40)
This relationship is satisfied if k satisfies
| k | | k | 2 .
(41)
That is, the assumption of slow-varying k implies that the WKB approximation is reasonable. This can also be recognized as the second derivative assumption [Eq. (30)] made above in the derivation of the zeroth-order WKB ODE.

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 long-wave region. Here, the wavenumber varies slowly in space, therefore, the left-hand 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 short-wave region. Here, the right-hand 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.

Further interrogation of the strong area assumption of Eq. (27) is in order as it also regards the spatial variation of model parameters. Under the first approximation for C0 [Eq. (32)], the assumption becomes
| S S | | k | .
(42)
This assumption appears to be challenged at the base, where | S / S | may be large as scala area varies approximately exponentially.60 

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 low-frequency stimuli, hence, box models without tapering reasonably satisfy the WKB approximation across most of space and frequency.

There are a number of methods for arriving at WKB approximate solutions for the 2-D model. The first solution presented in this tutorial is chosen because of its emphasis of the relationship between the 1-D and 2-D models. Under the assumptions of the model (see Sec. III), transmembrane pressure satisfies the 2-D Laplace equation such that
2 p x 2 + 2 p z 2 = 0.
(43)
One classical method for solving the Laplace equation is separation of variables, where it is assumed that the transmembrane pressure can be written as a product of a function of only x and a function of only z:
p ( x , z ) = X ( x ) Z ( z ) .
If separation of variables were satisfied, the solution would be a linear combination of eigenfunctions with eigenvalue k. These are of the forms
p k = ( A cosh k x + B cosh k x ) ( C e jkz + D e jkz ) ,
(44)
p k = ( A e jkx + B e jkx ) ( C cosh k z + D sinh k z ) .
(45)
We know that the x-dependence of the solution should have the form of a wave, therefore, Z should have the hyperbolic trigonometric form observed in Eq. (45):
Z ( z ) = C cosh k z + D sinh k z .
Plugging in the boundary condition at the outer wall [Eq. (7)] gives
Z ( h ) = k [ C sinh k h + D cosh k h ] = 0 ,
yielding the relationship C = D / tanh k h. A hyperbolic trigonometric identity results in
Z ( z ) = D sinh k h ( sinh k h sinh k x cosh k h cosh k x ) = D sinh k h cosh [ k ( z h ) ] .
(46)
Separation of variables has already been broken as h depends on x, but this nonetheless gives motivation for writing the form of the solution as
p ( x , z ) = cosh [ k ( z h ) ] X ( x ) ,
where k may also vary in x. In this form, the solution satisfies the boundary condition at the outer wall.
Here, I will make the first of two WKB approximations by assuming that (1) the form of X is that of Eq. (2) with δ = 1, and (2) the WKB assumption (parameters vary slowly in x relative to their own magnitudes) holds such that x-derivatives of cosh [ k ( z h ) ] are small. The Laplace equation to zeroth-order becomes
C 0 2 cosh [ k ( z h ) ] X ( x ) + k 2 cosh [ k ( z h ) ] X ( x ) = 0 ,
where the second term is the second z-derivative of p. Just as in the 1-D case, this gives
C 0 = ± j 0 x k ( ξ ) d ξ .
Assuming that the basal-traveling wave is negligible, the pressure is
p ( x , z ) = A ( x ) cosh [ k ( z h ) ] e j 0 x k ( ξ ) d ξ ,
(47)
where A(x) is some unknown function.
A dispersion relation is still needed so that k can be solved for. To do this, recall the boundary condition at z = 0 [Eq. (8)]. Pressure and velocity potential are related by p = 2 ρ j ω ϕ [Eq. (14)], therefore, their first derivatives in z give
w ̇ = 1 2 ρ j ω p z ,
or by plugging in Eq. (47),
w ̇ = 1 2 ρ j ω A ( x ) k sinh [ k ( z h ) ] e j 0 x k ( ξ ) d ξ .
(48)
At z = 0, the ratio of velocity and pressure is the negative of the OCC admittance, Y O C. That is,
Y O C = 1 2 j ω ρ A ( x ) k sinh [ k h ] e j 0 x k ( ξ ) d ξ A ( x ) cosh [ k h ] e j 0 x k ( ξ ) d ξ = k tanh [ k h ] 2 j ω ρ ,
giving the dispersion relation
k tanh k h = 2 j ω ρ Y O C .
(49)
This dispersion relation61 is transcendental and does not possess a unique solution for k (the implications of this will be covered in detail in Sec. VII). Equation (49) is also independent of A, meaning that it will be valid for any approximation of p in the form of Eq. (47).
Solving Eq. (49) for impedance gives
Z O C = j ω 2 ρ k tanh k h .
This allows for an attractive interpretation of the impact of the impedance on the traveling wave. Due to the leading j ω, this appears similar to a mass. In particular, defining the effective height of the fluid as
h e ( k ) = 1 k tanh k h , Z O C = 2 j ω ρ h e ,
(50)
the impedance at the OCC is that of a column of fluid with this effective height. It should be noted that while this analogy is useful, this “mass” and “height” are generally complex valued and only approximately real at lower frequencies (see discussion of the long-wave solution in Sec. V C).
To determine A(x) in Eq. (47), recall that the average pressure in the 2-D model must satisfy the Webster horn equation [Eq. (22)]. The average pressure is
p ¯ ( x ) = 1 h ( x ) 0 h ( x ) A ( x ) cosh [ k ( z h ) ] e j 0 x k ( ξ ) d ξ d z = 1 k ( x ) h ( x ) A ( x ) sinh [ k ( x ) h ( x ) ] e j 0 x k ( ξ ) d ξ .
(51)
The pressure focusing factor, α = p ( x , 0 ) / p ¯ ( x ), is required to find k 2 D in Eq. (23) and can now be found using Eq. (47):
α ( x ) = k ( x ) h ( x ) tanh [ k ( x ) h ( x ) ] .
(52)
This is independent of A, which means that it will be valid for any p in the form of Eq. (47).
Plugging this into Eq. (23), k 2 D 2 can be found to be
k 2 D 2 = 2 j ω Y O C ρ k h h tanh [ k h ] ,
but by Eq. (49), this simplifies directly to
k 2 D 2 = k 2 .
In Sec. IV, zeroth- and first-order WKB approximations for solutions to the Webster horn equation were derived. This gives an approximate formula for average pressure by directly copying Eq. (38):
p ¯ ( x ) = P O W S 0 k 0 S ( x ) k ( x ) e j 0 x k ( ξ ) d ξ .
(53)
Equating this with the earlier expression for p ¯ in Eq. (51), it is possible to solve for A(x):
A ( x ) = P O W k ( x ) h ( x ) sinh [ k ( x ) h ( x ) ] S 0 k 0 S ( x ) k ( x ) .
(54)
Finally, after this second application of a WKB approximation, a 2-D equation for pressure has been derived:
p ( x , z ) = P O W k ( x ) h ( x ) sinh [ k ( x ) h ( x ) ] S 0 k 0 S ( x ) k ( x ) cosh [ k ( x ) ( z h ( x ) ) ] e j 0 x k ( ξ ) d ξ .
(55)
In conjunction with the dispersion relation of Eq. (49), this allows solution for pressure or velocity throughout the scala.

The above derivation arrives at Eq. (55) through two consecutive applications of the WKB approximation and neatly piggy-backs off of 1-D results for average pressure. However, this formula is not the only solution referred to in the literature as “the WKB solution” for a 2-D model.

Various alternate approximation methods arrive at the following equation for pressure in a box model:
p ( x , z ) = P O W k 0 h cosh [ k ( x ) h ] tanh [ k 0 h ] tanh [ k 0 h ] + k 0 h sech 2 [ k 0 h ] tanh [ k ( x ) h ] + k ( x ) h sech 2 [ k ( x ) h ] × cosh [ k ( x ) ( z h ) ] e j 0 x k ( ξ ) d ξ
(56)
(e.g., Refs. 21, 27, and 28). It is clear that this has the form of Eq. (47), which means that this solution shares the same effective height, dispersion relation, and pressure focusing factor as the solution derived above [Eqs. (49), (50), and (52)].

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 Euler-Lagrange equations in a lossless box model (i.e., ZOC 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 lower-order approximation of Eq. 55.29–31,41,63 Differences in the behaviors between these solutions will be discussed in Sec. VIII.

It is also instructive to consider the behavior of the solution in the long-wave (small k, basal to best place/lower frequency than best frequency) and short-wave (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 : tanh a a if a is small and tanh a 1 if a is large.

The dispersion relations [Eq. (49)] in the long-wave and short-wave limits are, respectively,
k l w 2 = 2 j ω ρ Z O C h ,
(57)
k s w = 2 j ω ρ Y O C .
(58)
These are explicit solutions for the wavenumber in these regions, simplifying computation. Notably, klw is precisely the wavenumber from the 1-D Webster horn equation [Eq. (20)].
Considering the same limiting behavior for the pressure focusing factor [Eq. (52)] gives
α l w = 1 ,
(59)
α s w = k h .
(60)
This reinforces the realization that the long-wave approximation and 1-D solution are equivalent at z = 0. The effective height from Eq. (50) also has corresponding long- and short-wave approximations such that
h e , l w = 1 h k 2 ,
(61)
h e , s w = 1 k .
(62)

To visualize the differences between the long-wave, short-wave, and WKB approximations, one can observe the effective height as k varies. Figure 2 shows he for the three solutions across various values of positive real k with h = 1 mm. It can be observed that the WKB solution for he exhibits a continuous switch-off between the long- and short-wave approximations near the point where these solutions intersect. Behaviors of long- and short-wave velocity responses are discussed in Sec. VIII.

FIG. 2.

(Color online) Effective height as a function of real k for the long-wave and short-wave approximations, alongside the 2-D WKB solution. The scala height, h, is set to 1 mm.

FIG. 2.

(Color online) Effective height as a function of real k for the long-wave and short-wave approximations, alongside the 2-D WKB solution. The scala height, h, is set to 1 mm.

Close modal

Sections IV and V described derivations of explicit equations for pressure in 1-D or 2-D tapered box models via the WKB approximation. These formulas are valid where the model parameters do not change rapidly relative to their magnitudes and basal-traveling 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 1-D or 2-D cochlear model that describes the pressure at the OCC in response to a single-frequency stimulus. This is an infinite-dimensional object, living in the Hilbert space, H, of smooth functions mapping from I = [ 0 , L ] (the interval of along which the OCC is modeled to span) to .

Whereas the set of exact solutions is a subset of this infinite-dimensional space, they will likely have qualitatively similar traveling wave forms for various perturbations to parameters and boundary conditions. Thus, they are likely to be well-approximated as living in a lower-dimensional subspace containing functions that resemble cochlear traveling waves. This motivates the concept of a traveling wave subspace.

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 1-D 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 1-D box model WKB basis functions are the most commonly observed in the literature.10–15,44

By Eq. (38) with constant cross-sectional area, the apical-traveling wave is proportional to
W + = 1 k e j 0 x k ( ξ ) d ξ .
(63)
The basal-traveling wave has been ignored thus far in this tutorial. However, Eq. (37) implies that the form of the WKB approximate solution for the basal-traveling wave would differ from W + only by the sign in the exponential. I define
W = 1 k e j 0 x k ( ξ ) d ξ , x I .
(64)
The set β = { W + , W } H is linearly independent, and its span, W = span ( β ), is a 2-D subspace of H, which I will refer to as the WKB wave-space. Any function, p W, can be written in terms of its apical-traveling ( p +) and basal-traveling ( p ) components as
p ( x ) = p + ( x ) + p ( x ) = ψ + W + ( x ) + ψ W ( x ) , x I ,
(65)
where the coefficients ψ ± are complex-valued constants.
One can form a system of two equations in two variables by differentiating either side such that
p x = ψ + W + x + ψ W x .
(66)
Solution for the coefficients is neatly written in terms of the Wronskian of β,10,12,66 which is
D = det ( W + W W + W ) = 2 j .
(67)
With the Wronskian, the projections onto each basis function can be written as
p ± = P ± [ p ] = ψ ± W ± = ± W ± D ( W x W x ) p = 1 2 ( 1 ± j k 2 k 2 ± j k x ) p ,
(68)
where P ± represents the operators projecting functions in H onto W ±.67,68

Of course, any exact solution to the BVP will not live in W, hence, the values for ψ ± found through this formula will not be constant. Thus, the projections are merely approximations that are best if the derivatives of ψ ± are sufficiently small.69 

Having these projection operators, one can formulate a numerical method for determining the apical- and basal-traveling 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 short-wave solutions, long-wave solutions, or WKB solutions in a tapered box model.

One natural application of the projection described above is the study of reflections in the cochlea. The basal (+) and apical (–) reflection coefficients can be defined as
R ± ( x ) = p ( x ) p ± ( x ) .
(69)
In a model of the cochlea in which fluid pressure is driven at the stapes, a “perfectly efficient” cochlea would reflect no energy in the basal direction (this is assumed in the derivations of Secs. IV and V) and R + would be zero. In a cochlear model that exhibits some inefficiency, this will be a spatially varying complex-valued quantity. Some reasonable sources of such reflections include roughness in the OCC impedance or nonlinearity.

Conversely, one can consider how basal-traveling waves reflect toward the apex via R . 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 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 R +) such that basal-traveling 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 model10—an important result to keep in mind for the modeling of OAEs.

The WKB basis functions may also be used as an analytic tool in finding approximate solutions to related cochlear models. Once again, I will specify to 1-D box models with constant area. Starting with Eq. (19), the dynamics are governed by a wave equation with spatially varying wavenumber,
d 2 p d x 2 + k 2 ( x ) p = 0 ,
(70)
where I have replaced the partial derivatives with ordinary derivatives as time-dependence is implicit. This is a second-order linear homogeneous ODE for which W ± are approximate, linearly independent solutions. If they were truly solutions, β = { W + , W } would form a fundamental set for this ODE, and its general solution would be given by Eq. (65).
A corresponding nonhomogeneous ODE to Eq. (70) is given by
d 2 p d x 2 + k 2 ( x ) p = g ( x ) ,
(71)
where g is some nonzero function defined on I, which is known as the forcing function. Equation (70) is the complementary equation (or associated homogeneous equation) for this nonhomogeneous ODE. The theory of linear ODEs tells us that the general solution to Eq. (71), pgen, can be written as the sum of the general solution to the complementary equation, pc, and any particular solution to Eq. (71), pp.70 That is,
p gen = p c + p p .
The complementary solution, pc, is a linear combination of the functions in the complementary equation's fundamental set, which can be approximated by the set of WKB solutions, β. That is,
p c a + W + + a W , a ± .
(72)
The theory of variation of parameters70 then gives a particular solution in terms of the functions in the fundamental set and the forcing function,
p p = W D 0 x W + ( ξ ) g ( ξ ) d ξ W + D 0 x W ( ξ ) g ( ξ ) d ξ ,
where D is, once again, the Wronskian of the WKB functions that was already found to be 2j in Eq. (67).
This results in an approximate closed-form general solution for Eq. (71):
p gen = [ a + 1 2 j 0 x W ( ξ ) g ( ξ ) d ξ ] W + + [ a + 1 2 j 0 x W + ( ξ ) g ( ξ ) d ξ ] W = p + W + + p W .
(73)
The values of a ± are found through the boundary conditions. In particular, if a known pressure is applied at the stapes (x = 0), creating an initial apical-traveling wave, we would have a + = p 0 (known constant) and a  = 0. This gives the parameter-free closed-form solution for the stapes-driven nonhomogeneous 1-D model as
p = [ p 0 1 2 j 0 x W ( ξ ) g ( ξ ) d ξ ] W + + [ 1 2 j 0 x W + ( ξ ) g ( ξ ) d ξ ] W .
(74)

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 apical-traveling component is modulated by the basal-traveling wave weighted by g and vice versa.

In Sec. VI B, I discussed the application of projection onto WKB waves to approximating local reflection phenomena [Eq. (69)]. The application is natural in this analytic treatment as well given the p ± values in Eq. (73).

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.

Much work has been performed regarding the study of the impact of roughness on the impedance in generating intracochlear reflections.11–15,44 That is, if the smooth impedance, Zs, were modified by a small longitudinally varying perturbation,
Z ( x ) = Z s ( x ) + δ Z ( x ) ,
this would impact the wavenumber of the traveling waves in both directions [Eq. (20) in 1-D and Eq. (49) in 2-D]. One could also model this as a roughening of the wavenumber, where the smooth wavenumber would be ks and the roughened (squared) wavenumber would be
k 2 ( x ) = k s 2 ( x ) + δ k 2 ( x ) .
For example, δ k 2 ( x ) may be modeled as samples from independent identically distributed normal distributions at each x. The roughness could also be designed to depend on stimulus frequency, but this dependence will be left implicit as it will not impact the derivations.
Rewriting the wave equation in terms of the roughened wavenumber, we have
d 2 p d x 2 + [ k s 2 ( x ) + δ k 2 ( x ) ] p = 0 ,
which is, in fact, homogeneous and linear. However, δ k 2 is not necessarily differentiable—in fact, it ought not be because “rough” implies non-smooth. This precludes use of the WKB approximation in its current form as the WKB assumption [Eq. (41)] is not well posed.
Moving the δ k 2 term to the opposite side gives
d 2 p r d x 2 + k s 2 ( x ) p = δ k 2 p ,
which is still homogeneous as the right-hand side is proportional to the dependent variable, p. The strategy is to approximate the right-hand side as a p-independent forcing function. If δ k 2 is small, we can consider this right-hand term as a perturbation to the otherwise smooth, complementary response, pc, of Eq. (72). In the case that an apical-traveling wave of magnitude p0 is induced at the stapes, this results in the approximation
δ k 2 p δ k 2 p 0 W + .
This is a known p-independent function, allowing the ODE to be interpreted as approximately nonhomogeneous.
That is, it is in the form of Eq. (71) with g = δ k 2 p 0 W +. The roughened pressure solution can be given by substituting this forcing function for g in Eq. (74):
p = p 0 [ 1 + 1 2 j 0 x δ k 2 ( ξ ) W + ( ξ ) W ( ξ ) d ξ ] W + p 0 [ 1 2 j 0 x δ k 2 ( ξ ) W + 2 ( ξ ) d ξ ] W .
(75)
This solution facilitates computation of the reflection coefficients through Eq. (69) in terms of the roughness function.71 

In Secs. II–VI, I have described theoretical underpinnings for WKB solutions to 1-D and 2-D box and tapered box models. In this section, I discuss the challenges involved in implementation of the derived model equations in software.

The 1-D model poses no such difficulty as the WKB pressure equation and dispersion relation are explicit and in terms of elementary functions, but the dispersion relation of Eq. (49) presents a challenge in the 2-D case. This relation is transcendental and, generally, has infinitely many solutions for k in the complex plane. To standardize the language, the solution for k is reframed as a root-finding problem for the function
f ( z ) = z tanh z C ,
(76)
where
z = k h , C = 2 ρ hj ω Y O C .
(77)

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 jkx, the solution should possess a positive real part to correspond to an apical-traveling 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.

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 (x-resolution is 7 μ 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.

FIG. 3.

(Color online) Root loci for f(z) in response to four stimulus frequencies using parameters from Steele and Taber (Ref. 27; see Table I). For a single color, the most basal position corresponds to the smallest real root. As x increases in 7 μm increments, the root traverses a clockwise arc in the fourth quadrant, eventually arriving at the negative imaginary axis.

FIG. 3.

(Color online) Root loci for f(z) in response to four stimulus frequencies using parameters from Steele and Taber (Ref. 27; see Table I). For a single color, the most basal position corresponds to the smallest real root. As x increases in 7 μm increments, the root traverses a clockwise arc in the fourth quadrant, eventually arriving at the negative imaginary axis.

Close modal

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 ZOC 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 

TABLE I.

Parameters used in all simulations. Physical parameters are taken from Steele and Taber (Ref. 27).

Parameter Symbol Value
Mass  m  1.5 × 10 3 g/mm2 
Resistance  r  2 × 10 6 Ns/mm3 
Stiffness  s(x 10 e 0.2 x N/mm3 
Scala height  h  1 mm 
Cochlea length  L  35 mm 
Fluid density  ρ  10 3 g/mm3 
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 x-space  N  1024 
Points in z-space (finite difference method)  N/A  16 
Parameter Symbol Value
Mass  m  1.5 × 10 3 g/mm2 
Resistance  r  2 × 10 6 Ns/mm3 
Stiffness  s(x 10 e 0.2 x N/mm3 
Scala height  h  1 mm 
Cochlea length  L  35 mm 
Fluid density  ρ  10 3 g/mm3 
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 x-space  N  1024 
Points in z-space (finite difference method)  N/A  16 

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 root-finding 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 root-finding 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 N-length vector with resolution Δ x. I will write the estimate for k at position xn as k ̂ n , n = 1 , 2 , , 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 long-wave region. This motivates the initial approximation of k ̂ 1 = k l w ( x 1 ). This can be used as the initial value in a standard root-finding algorithm such as Newton-Raphson or the Muller method, which are likely to converge to the correct root.

Stepping further along in x, the long-wave 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 k n + 1. The intuitive estimate is to use the solution at xn, k ̂ n, as the starting point for the root-finding method at x n + 1 to find k ̂ n + 1.

Pseudocode for this algorithm using the Newton-Raphson method in z to compute the wavenumber is presented in Algorithm 1. The Newton-Raphson method requires the derivative of f, which is given by
f ( z ) = tanh z + z sech 2 z .
(78)
Recall that z = kh.

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

ALGORITHM 1.

Naive longitudinally stepping root-finding algorithm to determine the wavenumber at N different x positions using the Newton-Raphson method.

k ̂ i c k l w ( x 1 )     Initialize using long-wave k 

for n = 1 N do z h k ̂ i c

 
N is the number of steps in x space 

for m = 1 M do z z f ( z ; x n ) f ( z ; x n )

 
M is the number of Newton-Raphson iterations 
end for   
k ̂ n z / h k ̂ i c k n  Initial value for next step is current guess for k 
end for 
k ̂ i c k l w ( x 1 )     Initialize using long-wave k 

for n = 1 N do z h k ̂ i c

 
N is the number of steps in x space 

for m = 1 M do z z f ( z ; x n ) f ( z ; x n )

 
M is the number of Newton-Raphson iterations 
end for   
k ̂ n z / h k ̂ i c k n  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 switch-off between dominance of two modes22–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.

To account for the issues found in tracing the continuous root locus near the resonant position, one can introduce a discontinuity into the x-stepping algorithm. This is performed by changing the initial value for the root-finding 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.

To begin, the term C in the root-finding problem [see Eq. (77)] is approximated near the resonant position. Where the stiffness and mass cancel, the admittance is Y O C 1 / R d, a real conductance (where the d subscript denotes evaluation near the resonant position), i.e., C is purely negative imaginary. I define a new term γ such that
γ = R d 2 ρ h ω .
(79)
The new problem to solve becomes
z tanh z = j γ .
(80)
The solution to this transcendental equation can be approximated using the assumption that z is small in magnitude and lives close to the imaginary axis. This ensures that the chosen value of k will correspond to the dominant mode, falling off less rapidly than what would be found by tracing the continuous root locus. A derivation of this solution,76 relying on Taylor expansions, is given in brief by Viergever.21 It yields
k d π 2 h γ j π 2 h ( 1 γ 2 ) , γ = R d 2 ρ h ω .
(81)

The discontinuous x-stepping algorithm traces the continuous root locus up to some xd at which it is determined that the WKB assumption is being violated. This can be determined before simulation21 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 x-stepping method is paused once | k ̂ n k ̂ n 1 | / Δ x > T.

After this point, kd of Eq. (81) is used as an initial step in the root-finding algorithm. If | k ̂ n k ̂ n 1 | / Δ 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 1). At each subsequent step, it is determined whether k ̂ satisfies this threshold—it will eventually do so, at which point we continue the locus-tracing process along this second locus. Pseudocode for this algorithm is presented in Algorithm 2.

ALGORITHM 2.

Longitudinally stepping root-finding algorithm to determine the wavenumber at N different x positions, accounting for the discontinuity in the wavenumber.

k ̂ i c k l w ( x 1 )       Initialize using long-wave k 
for n = 1 N do      N is the number of steps in x space 
   z h k ̂ i c   

for m = 1 M do

   z z f ( z ) f ( z )

 
M is the number of Newton-Raphson iterations 
  end for   
  if | z / h k i c | Δ x < T then  Check for validity of WKB condition 
    k ̂ n z / h      Initial value for next step is current guess for k 
    k ̂ i c k n 
  else   
    k ̂ n NaN  Pressure and velocity should not be computed here 
    k ̂ i c k d  Guess for k after the discontinuity 
  end if   
end for   
k ̂ i c k l w ( x 1 )       Initialize using long-wave k 
for n = 1 N do      N is the number of steps in x space 
   z h k ̂ i c   

for m = 1 M do

   z z f ( z ) f ( z )

 
M is the number of Newton-Raphson iterations 
  end for   
  if | z / h k i c | Δ x < T then  Check for validity of WKB condition 
    k ̂ n z / h      Initial value for next step is current guess for k 
    k ̂ i c k n 
  else   
    k ̂ n NaN  Pressure and velocity should not be computed here 
    k ̂ i c k d  Guess for k after the discontinuity 
  end if   
end for   

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 long-wave approximation at every frequency-location pair, k ̂ = k l w. Then, k ̂ is plugged into Eq. (52) to find an the approximate pressure focusing factor, α ̂, and then α ̂ is plugged into Eq. (23) to find a new wavenumber approximation, k ̂. This is repeated for some number of iterations. Pseudocode for this algorithm is shown in Algorithm 3.

ALGORITHM 3.

Fixed point algorithm that updates a vector of k ̂ approximations by iteratively applying two relationships.

k ̂ k l w  Here, k ̂ is a vector with an index for each position 
for m = 1 M do  M is the number of fixed point iterations 
   α ̂ h k ̂ tanh k h  Pressure focusing vector update [Eq. (52)
   k ̂ 2 j ω ρ α ̂ h Z O C  Wavenumber update [Eq. (23)
  if R [ k ̂ ] < 0 then  Ensure the root is for an apical-traveling wave ( R gives the real part) 
    k ̂ k ̂ 
end if   
end for   
k ̂ k l w  Here, k ̂ is a vector with an index for each position 
for m = 1 M do  M is the number of fixed point iterations 
   α ̂ h k ̂ tanh k h  Pressure focusing vector update [Eq. (52)
   k ̂ 2 j ω ρ α ̂ h Z O C  Wavenumber update [Eq. (23)
  if R [ k ̂ ] < 0 then  Ensure the root is for an apical-traveling wave ( R gives the real part) 
    k ̂ k ̂ 
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.

Having developed the WKB approximate solutions as well as methods by which to find the wavenumber, it is now possible to observe the behavior of the modeled solutions. WKB solutions are compared to numerical results computed using the finite difference method of Neely.19 Physical quantities used here are from the 2-D box model of Steele and Taber.27 These, along with parameters used in wavenumber-finding algorithms, are provided in Table I. Mass, resistance, and stiffness terms contribute to the impedance according to
Z O C ( x ) = j ω m + r + s ( x ) j ω ,
(82)
where x is in units of mm.

In Sec. IV, I derived WKB solutions to the 1-D 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.

FIG. 4.

(Color online) Comparison of numerical solutions to the 1-D box model with zeroth- and first-order WKB approximate solutions in response to a 2 kHz stimulus.

FIG. 4.

(Color online) Comparison of numerical solutions to the 1-D box model with zeroth- and first-order WKB approximate solutions in response to a 2 kHz stimulus.

Close modal

It can be observed that the zeroth-order approximation overestimates the magnitude of the response near the peak and exhibits more phase accumulation than the numerical solution. On the other hand, the first-order 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.

To contextualize findings for the 2-D WKB solutions, it is useful to observe the performance of the long- and short-wave approximate solutions to the 2-D 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 non-monotonically 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 long-wave and short-wave solutions to the 2-D box model alongside a numerical solution. The long-wave 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.

FIG. 5.

(Color online) Comparison of numerical solutions to the 2-D box model with long- and short-wave approximate solutions in response to a 2 kHz stimulus.

FIG. 5.

(Color online) Comparison of numerical solutions to the 2-D box model with long- and short-wave approximate solutions in response to a 2 kHz stimulus.

Close modal

The short-wave 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 short-wave solution exhibiting a sign change in its group velocity just past the peak, is one dramatic consequence of this breakdown.

Before comparing 2-D WKB approximations to numerical results, it is first important to assess the methods for determining the wavenumber k in the 2-D case. This is performed by observing velocity responses at the OCC derived from the 2-D WKB approximation of Eq. (56), using three methods for finding the wavenumber: (1) Algorithm 1, an x-stepping algorithm that does not account for the discontinuity, amounting to following a root locus as in Fig. 3; (2) Algorithm 2, an x-stepping 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 x-stepping 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.

FIG. 6.

(Color online) Simulated velocity of the OCC in response to a 5.5 kHz stimulus using the parameters of Table I. Left-hand panels show the magnitude and phase responses for velocity computed using the 2-D WKB solutions with the dispersion relation solved via x-stepping methods. The method labeled as “discontinuous” accounts for the resonance by stopping computations when the derivative exceeds a threshold (Algorithm 2), and the method labeled as “continuous” simply follows a continuous root locus (Algorithm 1). The followed loci are displayed in the top-right panel, where the dashed arrow indicates the jump in the discontinuous method when entering the cutoff region. The bottom-right panel shows f(kh), which should be identically zero at roots.

FIG. 6.

(Color online) Simulated velocity of the OCC in response to a 5.5 kHz stimulus using the parameters of Table I. Left-hand panels show the magnitude and phase responses for velocity computed using the 2-D WKB solutions with the dispersion relation solved via x-stepping methods. The method labeled as “discontinuous” accounts for the resonance by stopping computations when the derivative exceeds a threshold (Algorithm 2), and the method labeled as “continuous” simply follows a continuous root locus (Algorithm 1). The followed loci are displayed in the top-right panel, where the dashed arrow indicates the jump in the discontinuous method when entering the cutoff region. The bottom-right panel shows f(kh), which should be identically zero at roots.

Close modal

Recall that these algorithms are designed to solve a root-finding problem for function f of Eq. (76). The root loci for f in the upper-right 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 bottom-right 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 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 x-stepping solution, but the phase accumulates more cycles.

FIG. 7.

(Color online) Same as that in Fig. 6 but also including the results of the fixed point algorithm as yellow curves in each subplot. Points to note about the fixed point algorithm include the erratic behavior of the root loci and failure to converge to a proper root near the transition between the peak and cutoff regions. Black arrows in the velocity response plots indicate positions at which this erratic behavior occurs.

FIG. 7.

(Color online) Same as that in Fig. 6 but also including the results of the fixed point algorithm as yellow curves in each subplot. Points to note about the fixed point algorithm include the erratic behavior of the root loci and failure to converge to a proper root near the transition between the peak and cutoff regions. Black arrows in the velocity response plots indicate positions at which this erratic behavior occurs.

Close modal

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 x-stepping 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 bottom-right 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 long-wave and short-wave 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

In Sec. VIII C, it was shown that the wavenumber-finding 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 lower-order solution of Eq. (55) and the higher-order solution of Eq. (56). Figure 8 shows solutions according to both of these equations alongside numerical solutions to the 2-D BVP.

FIG. 8.

(Color online) Comparison of numerical solutions to the 2-D box model with lower-order [Eq. (55)] and higher-order [Eq. (56)] WKB approximate solutions in response to a 2 kHz stimulus.

FIG. 8.

(Color online) Comparison of numerical solutions to the 2-D box model with lower-order [Eq. (55)] and higher-order [Eq. (56)] WKB approximate solutions in response to a 2 kHz stimulus.

Close modal

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 root-finding methods, as observed in Fig. 7.

The lower-order solution slightly overestimates tuning at the peak as compared to the higher-order 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.

The WKB approximation provides compact closed-form approximate solutions for 1-D and 2-D 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 apical-traveling wave. Through the method of WKB wave-space projection, the solutions also facilitate further interpretation of modeled responses as a superposition of apical- and basal-traveling 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 2-D 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 uniform-height 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 cell-Deiters 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 level-dependent 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.

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 DC020621-02).

The author has no conflicts to disclose.

No data were used in this tutorial article. Scripts used in generating the presented figures are available from the author upon reasonable request.

In this appendix, I provide a derivation of the higher-order “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 WKB-adjacent ansatz (but not precisely the WKB method).

The method consists of the following steps:

  1. 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

  2. write a form for the solution to this new PDE as
    A ( x , ζ ) cosh [ κ ( x ) ( H ζ ) ] 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.

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

  4. 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.

TABLE II.

Symbols introduced in the derivation of the model equations in the series solution approach along with their significance and units.

Symbol Significance Units
Z0  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 x-dependence of the PDE  Unitless 
ζ  Kz, Nondimensionalized transverse coordinate  Unitless 
H  Kh, Nondimensionalized scala height  Unitless 
Q ( x , ζ )  Pressure written in terms of the nondimensionalized transverse coordinate  Pa 
A ( x , ζ )  Auxiliary pressure variable that controls the magnitude of pressure at the OCC, to be solved for in the simplified BVP  Pa 
κ ( x )  Controls the x-dependence 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
Z0  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 x-dependence of the PDE  Unitless 
ζ  Kz, Nondimensionalized transverse coordinate  Unitless 
H  Kh, Nondimensionalized scala height  Unitless 
Q ( x , ζ )  Pressure written in terms of the nondimensionalized transverse coordinate  Pa 
A ( x , ζ )  Auxiliary pressure variable that controls the magnitude of pressure at the OCC, to be solved for in the simplified BVP  Pa 
κ ( x )  Controls the x-dependence 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 
Viergever begins with the 2-D box model Laplace equation BVP in P(x,z), and then performs a change of variables. To start, we define a reference impedance, Z0, which is some arbitrary constant. We also define f 2 ( x ) = Z 0 / Z O C ( x ) and a reference wavenumber, K 2 = 2 j ω ρ / h Z 0. Recalling that P = ρ ϕ ̇, the boundary condition at the OCC is
P z + h K 2 f 2 ( x ) P = 0 , z = 0.
(A1)
Note that K is not a function of x.
Further reparameterizing the z coordinate and defining a reparameterized pressure, Q, as
ζ = K z , H = K h , Q ( x , ζ ) = P ( x , z ) ,
(A2)
the BVP in terms of Q is
2 Q x 2 + K 2 2 Q ζ 2 = 0 ,
(A3)
Q ζ | ζ = H = 0 ,
(A4)
Q ζ | ζ = 0 + H f 2 ( x ) Q ( x , 0 ) = 0.
(A5)
The solution to the above PDE is artificially represented in a form resembling what the solutions are expected to be by intuition about the solutions of the Laplace equation. In particular, Q is written as
Q ( x , ζ ) = A ( x , ζ ; K ) e jKg ( x ) cosh [ κ ( x ) ( H ζ ) ] .
(A6)

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)].

Plugging this form of Q into the BVP, we can find an equivalent BVP in terms of A. Writing a ( x , ζ ) = κ ( x ) ( H ζ ) to simplify notation, we arrive at the following PDE and boundary conditions:
K 2 [ ( κ 2 g 2 ) A cosh a + 2 A ζ 2 cosh a 2 κ A ζ sinh a ] + j K [ g A cosh a + 2 g A cosh a x ] + 2 A cosh a x 2 = 0 ,
(A7)
A ζ | ζ = 0 κ A ( x , 0 ) tanh a ( x , 0 ) + H f 2 A ( x , 0 ) = 0 ,
(A8)
A ζ | ζ = H = 0.
(A9)

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
A formal power series solution in the form
A ( x , ζ ; K ) = A 0 ( x ) + n = 1 1 ( j K ) n A n ( x , ζ )
(A10)
is assumed with monotonic decrease in magnitude of terms and their derivatives in increasing n, and allows for termwise differentiation. This form of the solution is not quite the WKB ansatz, but it is the logarithm of such a solution with δ = j K. This is motivated by Keller's approach to surface waves.62 
This ansatz is plugged into the PDE for A in Eq. (A7), resulting in a system of infinitely many PDEs of which we consider only the PDEs that include A0 and A1 (justified by the assumption that the terms and their derivatives decrease monotonically). The resulting system of differential equations is
g 2 ( x ) = κ 2 ( x ) ,
(A11)
cosh a 2 A 1 ζ 2 2 κ sinh a A 1 ζ = g A 0 cosh a + 2 g A 0 cosh a x .
(A12)
Application of boundary conditions results in
A 1 ζ | ζ = H = 0 ,
(A13)
A 1 ζ | ζ = 0 = 0 ,
(A14)
κ tanh κ H = H f 2 .
(A15)
Equation (A15) resembles the dispersion relation derived from the WKB method in Sec. V.
Equation (A11) is solved by
g ( x ) = ± 0 x κ ( ξ ) d ξ + C
(A16)
for arbitrary constant C. This resembles the characteristic WKB phase term.
3. Finding a first approximation for pressure
Solving Eq (A12) is nontrivial as it contains A0 and A1. Solution for A0 requires clever substitutions.78 I find
A 0 = C ( κ H + sinh κ H cosh κ H ) 1 / 2
(A17)
for arbitrary C.

Theoretically, this facilitates solution for An 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).

Ignoring An for n 1 gives a first approximation for Q by putting A A 0. Using Eq. (A16) for g, there are two possible solutions:
Q ± ( x , ζ ) = C ± ( κ H + sinh κ H cosh κ H ) 1 / 2 e ± j K 0 x κ ( ξ ) d ξ cosh [ κ ( x ) ( H ζ ) ] + O ( 1 / K ) .
(A18)

The reference constant, K, was defined as K 2 = 2 j ρ ω / h Z 0, where Z0 was a second reference constant such that f 2 = Z 0 Y O C. Because Z0 was entirely arbitrary, I am free to choose Z 0 = 2 j ρ ω h 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).

I define k = 1 mm 1 × κ. In this light, the κ relation from Eq. (A15) becomes
k tanh k h = 2 j ρ ω Y O C .
(A19)
This is precisely the dispersion relation derived through the WKB method [Eq. (49)], where k is the wavenumber with units mm−1.
The first approximation for pressure with arbitrary constants, C ±, is now
P ( x , z ) = ( k h + sinh k h cosh k h ) 1 / 2 cosh [ k ( x ) ( h z ) ] [ C + e j 0 x k ( ξ ) d ξ + C e j 0 x k ( ξ ) d ξ ] .
(A20)
To find the constants, the two x-boundary conditions are used:
1 h 0 h P ( 0 , z ) d z = P O W , P x | x = L = 0 ,
(A21)
where L is the length of the cochlea, and POW is the average pressure at the stapes.
After some computation,79 assuming that the backward traveling wave is negligible, we achieve
P ( x , z ) = P O W k 0 h cosh k ( x ) h tanh k 0 h k 0 h sech 2 k 0 h + tanh k 0 h k ( x ) h sech 2 k ( x ) h + tanh k ( x ) h cosh [ k ( x ) ( h z ) ] e j 0 x k ( ξ ) d ξ .
(A22)
This is precisely Eq. (56).
1.
R. L.
Wegel
and
C. E.
Lane
, “
The auditory masking of one pure tone by another and its probable relation to the dynamics of the inner ear
,”
Phys. Rev.
23
(
2
),
266
285
(
1924
).
2.
J.
Zwislocki
, “
Theory of cochlear mechanics
,”
Hear. Res.
2
(
3–4
),
171
182
(
1980
).
3.
L. C.
Peterson
and
B. P.
Bogert
, “
A dynamical theory of the cochlea
,”
J. Acoust. Soc. Am.
22
(
3
),
369
381
(
1950
).
4.
G.
Zweig
,
R.
Lipes
, and
J. R.
Pierce
, “
The cochlear compromise
,”
J. Acoust. Soc. Am.
59
(
4
),
975
982
(
1976
).
5.
G.
Zweig
, “
Finding the impedance of the organ of Corti
,”
J. Acoust. Soc. Am.
89
(
3
),
1229
1254
(
1991
).
6.
C. A.
Shera
and
G.
Zweig
, “
A symmetry suppresses the cochlear catastrophe
,”
J. Acoust. Soc. Am.
89
(
3
),
1276
1289
(
1991
).
7.
S. J.
Elliott
,
E. M.
Ku
, and
B.
Lineton
, “
A state space model for cochlear mechanics
,”
J. Acoust. Soc. Am.
122
(
5
),
2759
2771
(
2007
).
8.
M. A.
Viergever
, “
Asymmetry in reflection of cochlear waves
,”
Aud. Freq. Sel.
119
,
31
38
(
1986
).
9.
E.
de Boer
,
C.
Kaernbach
,
P.
König
, and
T.
Schillen
, “
Forward and reverse waves in the one-dimensional model of the cochlea
,”
Hear. Res.
23
(
1
),
1
7
(
1986
).
10.
C. A.
Shera
and
G.
Zweig
, “
Reflection of retrograde waves within the cochlea and at the stapes
,”
J. Acoust. Soc. Am.
89
(
3
),
1290
1305
(
1991
).
11.
C. L.
Talmadge
,
A.
Tubis
,
G. R.
Long
, and
P.
Piskorski
, “
Modeling otoacoustic emission and hearing threshold fine structures
,”
J. Acoust. Soc. Am.
104
(
3
),
1517
1543
(
1998
).
12.
C. L.
Talmadge
,
A.
Tubis
,
G. R.
Long
, and
C.
Tong
, “
Modeling the combined effects of basilar membrane nonlinearity and roughness on stimulus frequency otoacoustic emission fine structure
,”
J. Acoust. Soc. Am.
108
(
6
),
2911
2932
(
2000
).
13.
C. A.
Shera
,
A.
Tubis
, and
C. L.
Talmadge
, “
Coherent reflection in a two-dimensional cochlea: Short-wave versus long-wave scattering in the generation of reflection-source otoacoustic emissions
,”
J. Acoust. Soc. Am.
118
(
1
),
287
313
(
2005
).
14.
A.
Moleti
and
R.
Sisto
, “
Localization of the reflection sources of stimulus-frequency otoacoustic emissions
,”
J. Assoc. Res. Otolaryngol.
17
(
5
),
393
401
(
2016
).
15.
R.
Sisto
,
A.
Moleti
, and
C. A.
Shera
, “
On the spatial distribution of the reflection sources of different latency components of otoacoustic emissions
,”
J. Acoust. Soc. Am.
137
(
2
),
768
776
(
2015
).
16.
W. M.
Siebert
, “
Ranke revisited—A simple short-wave cochlear model
,”
J. Acoust. Soc. Am.
54
,
282
(
1973
).
17.
J. B.
Allen
, “
Two-dimensional cochlear fluid model: New results
,”
J. Acoust. Soc. Am.
61
(
1
),
110
119
(
1977
).
18.
J. B.
Allen
and
M. M.
Sondhi
, “
Cochlear macromechanics: Time domain solutions
,”
J. Acoust. Soc. Am.
66
(
1
),
123
132
(
1979
).
19.
S. T.
Neely
, “
Finite difference solution of a two-dimensional mathematical model of the cochlea
,”
J. Acoust. Soc. Am.
69
,
1386
1391
(
1981
).
20.
C. R.
Steele
and
C. E.
Miller
, “
An improved WKB calculation for a two-dimensional cochlear model
,”
J. Acoust. Soc. Am.
68
(
1
),
147
148
(
1980
).
21.
M. A.
Viergever
,
Mechanics of the Inner Ear: A Mathematical Approach
(
Technische Hogeschool
,
the Netherlands
,
1980
).
22.
L.
Watts
, “
Cochlear mechanics: Analysis and analog VLSI
,” Ph.D. thesis,
California Institute of Technology
,
Pasadena, CA
,
1993
.
23.
L.
Watts
, “
The mode-coupling Liouville–Green approximation for a two-dimensional cochlear model
,”
J. Acoust. Soc. Am.
108
(
5
),
2266
2271
(
2000
).
24.
S. J.
Elliott
,
G.
Ni
,
B. R.
Mace
, and
B.
Lineton
, “
A wave finite element analysis of the passive cochlea
,”
J. Acoust. Soc. Am.
133
(
3
),
1535
1545
(
2013
).
25.
J.
Lighthill
, “
Energy flow in the cochlea
,”
J. Fluid Mech.
106
,
149
213
(
1981
).
26.
J.
Lighthill
,
Waves in Fluids
(
Cambridge University Press
,
Cambridge, UK
,
1978
).
27.
C. R.
Steele
and
L. A.
Taber
, “
Comparison of WKB and finite difference calculations for a two-dimensional cochlear model
,”
J. Acoust. Soc. Am.
65
(
4
),
1001
1006
(
1979
).
28.
H.
Duifhuis
, “
Cochlear macromechanics
,” in
Mechanics of Hearing
(
Springer
,
New York
), Chap. 6, pp.
189
211
(
1988
).
29.
R.
Sisto
,
D.
Belardinelli
, and
A.
Moleti
, “
Fluid focusing and viscosity allow high gain and stability of the cochlear response
,”
J. Acoust. Soc. Am.
150
(
6
),
4283
4296
(
2021
).
30.
R.
Sisto
,
D.
Belardinelli
,
A.
Altoè
,
C. A.
Shera
, and
A.
Moleti
, “
Crucial 3-D viscous hydrodynamic contributions to the theoretical modeling of the cochlear response
,”
J. Acoust. Soc. Am.
153
(
1
),
77
86
(
2023
).
31.
C. A.
Shera
and
A.
Altoè
, “
Otoacoustic emissions reveal the micromechanical role of organ-of-Corti cytoarchitecture in cochlear amplification
,”
Proc. Natl. Acad. Sci. U.S.A.
120
(
41
),
e2305921120
(
2023
).
32.
E. S.
Olson
, “
Intracochlear pressure measurements related to cochlear tuning
,”
J. Acoust. Soc. Am.
110
(
1
),
349
367
(
2001
).
33.
This concept has been of particular interest in cochlear mechanics since Olson showed that pressure varies rapidly in space within the scala and is tuned near the OCC.
34.
J.
Liouville
, “
Extrait d'un mémoire sur le développement des fonctions en séries dont les différents termes sont assujettis satisfaire une méme équation différentielle linéaire, contenant un paramètre variable
” (“Extract from a dissertation on the development of functions in series whose different terms are subject to satisfying the same linear differential equation, containing a variable parameter”),
J. Math. Pures Appl.
2
,
16
35
(
1837
).
35.
G.
Green
, “
On the motion of waves in a variable canal of small depth and width
,”
Trans. Cambridge Philos. Soc.
6
,
457
462
(
1837
).
36.
J.
Mathews
and
R. L.
Walker
,
Mathematical Methods of Physics
(
Addison-Wesley
,
Boston, MA
,
1970
).
37.
R. B.
Dingle
,
Asymptotic Expansions Their Derivation and Interpretation
(
Academic
,
New York
,
1975
).
38.
S.
Winitzki
, “
Cosmological particle production and the precision of the WKB approximation
,”
Phys. Rev. D
72
(
10
),
104011
(
2005
).
39.
E.
de Boer
, “
Auditory physics. physical principles in hearing theory. II
,”
Phys. Rep.
105
(
3
),
141
226
(
1984
).
40.
E.
de Boer
and
M.
Viergever
, “
Wave propagation and dispersion in the cochlea
,”
Hear. Res.
13
(
2
),
101
112
(
1984
).
41.
A.
Altoè
,
J. B.
Dewey
,
K. K.
Charaziak
,
J. S.
Oghalai
, and
C. A.
Shera
, “
Overturning the mechanisms of cochlear amplification via area deformations of the organ of Corti
,”
J. Acoust. Soc. Am.
152
(
4
),
2227
2239
(
2022
).
42.
F.
Mammano
and
R.
Nobili
, “
Biophysics of the cochlea: Linear approximation
,”
J. Acoust. Soc. Am.
93
(
6
),
3320
3332
(
1993
).
43.
E.
De Boer
, “
Short waves in three-dimensional cochlea models: Solution for a ‘block’ model
,”
Hear. Res.
4
(
1
),
53
77
(
1981
).
44.
R.
Sisto
,
A.
Moleti
, and
C. A.
Shera
, “
Cochlear reflectivity in transmission-line models and otoacoustic emission characteristic time delays
,”
J. Acoust. Soc. Am.
122
(
6
),
3554
3561
(
2007
).
45.
Y.-J.
Yoon
,
C. R.
Steele
, and
S.
Puria
, “
Feed-forward and feed-backward amplification model from cochlear cytoarchitecture: An interspecies comparison
,”
Biophys. J.
100
,
1
10
(
2011
).
46.
M.
Robnik
and
V. G.
Romanovski
, “
Some properties of WKB series
,”
J. Phys. A: Math. Gen.
33
(
28
),
5093
5104
(
2000
).
47.
In reality, the series often diverges and terms will begin to increase after a certain order of m. This limits the precision of the approximation but will not come into play in this tutorial as I will never exceed m = 1. Detailed discussion is presented in Ref. 38.
48.
G.
Von Békésy
and
E. G.
Wever
,
Experiments in Hearing
(
McGraw-Hill
,
New York
,
1960
).
49.
Reissner's membrane, separating SV and scala media, is assumed to be mechanically invisible, therefore, we need not distinguish between all three chambers.
50.
This assumption can be relaxed and asymmetric scalae can be assumed by finding some effective area. Results are largely unchanged.
51.
M. A.
Viergever
, “
On the physical background of the point-impedance characterization of the basilar membrane in cochlear mechanics
,”
Acustica
39
,
292
297
(
1978
).
52.
L. J.
Kanis
and
E.
de Boer
, “
Self-suppression in a locally active nonlinear model of the cochlea: A quasilinear approach
,”
J. Acoust. Soc. Am.
94
(
6
),
3199
3206
(
1993
).
53.
L. J.
Kanis
and
E.
de Boer
, “
Comparing frequency-domain with time-domain solutions for a locally active nonlinear model of the cochlea
,”
J. Acoust. Soc. Am.
100
(
4
),
2543
2546
(
1996
).
54.
The focus of this tutorial is on linear models, but nonlinearity can be implemented in a number of ways (e.g., quasilinear method, formulation of a nonhomogeneous cochlear model, etc.).
55.
E.
Poisson
,
A.
Pound
, and
I.
Vega
, “
The motion of point particles in curved spacetime
,”
Living Rev. Relativ.
14
,
7
(
2011
).
56.
Retarded Green's functions are tools used in the study of general relativity that are often challenging to write in closed form. These are not to be confused with standard Green's function methods used in cochlear mechanics research, which generally give implicit solutions.
57.
T.
Zamorski
and
R.
Wyrzykowski
, “
Approximate methods for the solution of the equation of acoustic wave propagation in horns
,”
Arch. Acoust.
6
(
3
),
273
285
(
1981
).
58.
V.
Salmon
, “
Generalized plane wave horn theory
,”
J. Acoust. Soc. Am.
17
(
3
),
199
211
(
1946
).
59.
Briefer derivations of the WKB approximate solution to the horn equation appear in Refs. 57 and 58 but are generally less detailed than the presentation in this tutorial.
60.
W.
Plassmann
,
W.
Peetz
, and
M.
Schmidt
, “
The cochlea in gerbilline rodents
,”
Brain, Behav. Evol.
30
(
1-2
),
82
102
(
1987
).
61.
A note on terminology: This equation is often called the eikonal equation in literature due to an analogy to geometric optics, but this language is imprecise. Eikonal equations are a class of differential equations that appear elsewhere in the above derivations, and this algebraic expression is not such an equation. Instead, I will simply use the term dispersion relation as it is more descriptive and precise.
62.
J. B.
Keller
, “
Surface waves on water of non-uniform depth
,”
J. Fluid Mech.
4
(
6
),
607
614
(
1958
).
63.
A.
Altoè
and
C. A.
Shera
, “
The cochlear ear horn: Geometric origin of tonotopic variations in auditory signal processing
,”
Sci. Rep.
10
(
1
),
20528
(
2020
).
64.
O. F.
Ranke
, “
Theory of operation of the cochlea: A contribution to the hydrodynamics of the cochlea
,”
J. Acoust. Soc. Am.
22
(
6
),
772
777
(
1950
).
65.
Notably, long- and short-wave solutions are not particular to the WKB approximation—for historical long-wave solutions, see the early work of Peterson and Bogert (Ref. 3) and for historical short-wave solutions, see Refs. 16 and 64.
66.
It should be noted that this Wronskian (also called the Wronskian determinant in some works) will be distinct for different choices of scaling factors for the basis functions. The projections will be identical.
67.
Throughout this article, functions are written in terms of wavenumber k. Wavelength is used in other works (Ref. 10), causing a difference in the formulas presented here, resolved simply by the chain rule of differentiation and the relationship between wavenumber and wavelength.
68.
A more intricate treatment of these derivations can be found at https://github.com/brian-lance/wkb-derivations (Last viewed January 10, 2024).
69.
A metric has been developed by Mathews and Walker and repeated by Shera and Zweig that quantifies the error through the size of these derivatives (see Refs. 10 and 36). This is done by representing the basis functions as exact solutions to a similar differential equation, W + k 2 ( x ) ( 1 + ϵ ) W = 0, where a smaller ϵ represents better approximations to the actual BVP.
70.
D. G.
Zill
,
W. S.
Wright
, and
M. R.
Cullen
,
Differential Equations with Boundary-Value Problems
, 8th ed. (
Brooks-Cole
,
Belmont, MA
,
2013
).
71.
The formula is similar to Eq. (74) of Ref. 12, except for (1) a negative sign and (2) complementary integral bounds. These discrepancies cancel if the integral over the interval is zero.
72.
This one-mode assumption is worthy of some scrutiny, and the works of Watts (Refs. 22 and 23) and Elliott (Ref. 24) have covered multimode solutions and implications thereof.
73.
E.
de Boer
and
M.
Viergever
, “
Validity of the Liouville-Green (or WKB) method for cochlear mechanics
,”
Hear. Res.
8
(
2
),
131
155
(
1982
).
74.
Note on terminology: Some modelers have described this point as the critical layer, owing its name to a more general theory of critical layer absorption (Refs. 25 and 26).
75.
L.
Robles
and
M. A.
Ruggero
, “
Mechanics of the mammalian cochlea
,”
Physiol. Rev.
81
(
3
),
1305
1352
(
2001
).
76.
An intricate treatment can be found at https://github.com/brian-lance/wkb-derivations (Last viewed January 10, 2024).
77.
A.
Latif
,
Banach Contraction Principle and Its Generalizations
(
Springer International
,
Cham
, 2014), pp.
33
64
.
78.
This is outlined at https://github.com/brian-lance/wkb-derivations (Last viewed January 10, 2024).
79.
These can be found at https://github.com/brian-lance/wkb-derivations (Last viewed January 10, 2024).