Effective-Spring Model of Tympanic Response in Archosaurs

Abstract

Whereas for smaller animals the eardrums are well-characterized as excitable membranes or drums, some animals such as several archosaurs feature, as a first approximation, a rather stiff elastic shell supported by an elastic ring. Mathematically, the theory of plates and shells is applicable but its governing equations overly complicate the modeling. Here the notion of tympanic structure is introduced as a generalization of “ordinary” tympanic membranes so as to account for sound perception as it occurs in archosaurs, such as birds and crocodilians. A mathematical model for the tympanic structure in many archosaurs called two-spring model implements this notion. The model is exactly soluble and solutions are presented in closed form and as a series expansion. Special emphasis is put onto offering an easy-to-apply model for describing experiments and performing numerical studies. The analytic treatment is supplemented by a discussion of the applicability of the two-spring model in auditory research. An elasticity-theoretic perspective of the two-spring model is given in the Appendix.

Share and Cite:

Heider, D. and van Hemmen, J. (2019) Effective-Spring Model of Tympanic Response in Archosaurs. Open Journal of Biophysics, 9, 21-50. doi: 10.4236/ojbiphy.2019.91003.

1. Introduction

1.1. Internally Coupled ears (ICE)

More than half of the land-living vertebrates possess an air-filled cavity connecting left and right eardrums. That is, they possess internally coupled ears, for short ICE [1] . Figure 1 shows the evolution of hearing in vertebrates using ICE [2] [3] [4] . In the simplest models [5] [6] of this widespread mechanism the two eardrums are an elastic structure interfacing the external acoustic environment

Figure 1. Left: Evolution of vertebrate hearing systems and animals using ICE. Courtesy of Anupam P. Vedurmudi for the full graphic, adapted from van Hemmen et al. [1] . Right: The elastic shell of a crocodilian (youngster’s) eardrum. The black bar sets a length scale of 1 cm. The diameter and height of the shell structure are 8 mm and 2 mm, respectively. The extracolumellar lever is located in the upper right quadrant of the picture. Its geometry resembles a spherical sector. Picture courtesy of Bruce A. Young.

with the interaural cavity, functioning as an acoustic wave guide [7] [8] . The interaural cavity ensures that the membranes’ vibrations are not independent but mutually coupled. For neuronally less developed animals, the coupling provides an easy-to-realize mechanism to generate a more precise directional hearing than without ICE; for early attempts at an explanation, see Autrum [9] [10] .

Although the mechanism of internally coupled ears seems in principio universal for numerous animals, the biological realization differs from species to species and animal to animal [1] . Aquatic frogs such as Xenopus [11] employ plate-like eardrums bounding an air-filled cavity whereas several birds, in particular chicken, use a flexible membrane in a cavity filled with, surprisingly, not air [3] . Lizards [5] [6] [12] use flexible membranes connected by a large cylindrical cavity whereas crickets employ, possibly to comply with their exoskeletal anatomy, a complicated tracheal system featuring two membranes connected through a complicated cavity system. The latter is partitioned itself into two symmetric halves by another elastic structure called septum [13] .

An increase in the area covered by the eardrums leads to a lowering of the fundamental frequency in a “clamped membrane model” because the first eigenvalue of the Laplacian (with Dirichlet boundary conditions) scales as the inverse of the square of the area covered by the eardrums. Thus a small eardrum area seems favorable to these animals. Anatomically, these requirements are realized in some archosaurs through a self-supporting spherical shell supported by an elastic ring. In reptiles and birds, an elastic ring occurs far more often [14] in the middle ear than one might naively expect.

The construction is shown schematically in Figure 2. The spherical shell itself may vibrate due to bending and is shown in vitro in Figure 1. Its net displacement in vertical direction leads to an induced vibration of the ring supporting the spherical shell elastically. This means that the notion of a tympanic membrane naturally breaks down because the shell is far less flexible than the elastically

Figure 2. Left: Schematic physical decomposition of the vibration of the tympanic structure due to an incident pressure signal. Right: Conceptual representation of the two-spring model. The pressure signal causes a uniform displacement of the spherical shell or, equivalently, an elastic deformation xr of the supporting ring underneath. The vibration of the ring leads, once the shell is not a hemisphere, to a an excess torque applied to the boundary of the shell and thus to shell bending, xs, and has been incorporated through a coupling of the xr-oscillator to the xs-oscillator through a massless spring of spring constant k0.

supporting ring so that the latter is responsible most for the fundamental frequency of the “eardrum”. It needs to be constantly borne in mind that the amplitudes are in the nm range whereas the ring has a thickness of mm so that the difference is at least 5 orders of magnitude. Physical auditory research [13] currently lacks both adequate jargon and formalism to capture these phenomena and indicate the direction of useful modeling. This is the gap we aspire to fill with the present treatment.

1.2. Tympanic Structures

A tympanic structure is defined to be the (finite) set of elastic elements that 1) interface the interaural cavity of the animal under consideration with its external acoustic environment, 2) by virtue of mutual elastic coupling of the constituents function as a single elastic structure, and 3) respond in this form at least locally to auditory signals. We call a given tympanic structure N-constituent tympanic structure, if the set of elastic elements in the sense of the previous definition has cardinality N. An elastic element of a tympanic structure is also called a constituent.

2. Two-Spring Model

2.1. Heuristic Derivation of the Two-Spring Model

The two-spring model is based on two simplifying assumptions. First, we are only interested in the impact that the vibration of a tympanic structure has on the volume in the interaural cavity. Vibrations give rise to a local mass density change in the cavity and thus trigger the formation of a pressure wave traveling between the, usually two, tympanic structures [1] [5] [6] . We refer to Howe [8] [7] for an introductory treatment of sound generated by vibrating elastic structures, to Timoshenko [15] [16] for a detailed treatment of elasticity, and to van Hemmen and Leibold [17] for elasticity theory as applied to biological membranes.

Let ( q 1 , q 2 , q 3 ) be coordinates such that q 3 = 0 and ( q 1 , q 2 ) Γ define the equilibrium position of an elastic element where Γ 2 is bounded and suitably regular domain. (For two eardrums there is, say, another element at q 3 = L > 0 .) Furthermore, let u be the displacement, i.e., a small u Area ( Γ ) perturbation from the former equilibrium. Solving elasticity equations for two-dimensional elastic structures is equivalent to the inversion of operators t 2 + O where O typically is the Laplacian Δ , its square Δ 2 , the bi-harmonic operator Δ bi 2 or linear combinations thereof [17] . Let us replace the displacement of the constituents of a tympanic structure by its average

x ( t ) 1 Area ( Γ ) Γ d 2 ( q 1 , q 2 ) u ( t , q 1 , q 2 ) . (1)

This approximation is called piston approximation [5] [12] . Assuming physical behavior of the full displacement u depending on the temporal variable t as well as the spatial coordinates q 1 , q 2 , we can apply the mean value theorem of integration. It ensures the existence of a ( q 1 , q 2 ) Γ such that u ( t , q 1 , q 2 ) = x ( t ) .

We now perform a partial Taylor series expansion in the spatial variables

u ( t , q 1 , q 2 ) = x ( t ) + R ( t , q 1 , q 2 ) (2)

where R ( t , q 1 , q 2 ) is Lagrange’s remainder,

R ( t , q 1 , q 2 ) = i = 1 2 ( q i q i ) 0 1 d s i u ( t , { ( 1 s ) q j + s q j } j ) . (3)

The remainder is absolutely bounded from above by firstly applying Cauchy’s inequality to the Euclidean inner product. We note that

( q 1 q 1 , q 2 q 2 ) 2 diam ( Γ ) sup { ( q 1 , q 2 ) ( q 1 , q 2 ) 2 : ( q 1 , q 2 ) , ( q 1 , q 2 ) Γ }

and that the function inside the integral is smaller than or equal to its maximum modulus. The remainder is now bounded by

| R | q 1 , q 2 ( t ) diam ( Γ ) max ( q 1 , q 2 ) Γ ( u ) 2 ( t ) . (4)

In particular, the piston approximation is dynamically accurate, if

max ( q 1 , q 2 ) Γ ( u ) 2 ( t ) | x | ( t ) diam ( Γ ) . (5)

Our first assumption is that the above inequality is fulfilled.

Next we define x ˜ and ω 2 > 0 through averaging over ( O ) [ u ] ,

ω 2 x ˜ ( t ) 1 Area ( Γ ) Γ d 2 ( q 1 , q 2 ) ( O ) [ u ] ( t , q 1 , q 2 ) . (6)

Since O 0 , i.e., O is positive-definite for all practical purposes, we obtain a square of an effective average frequency and an average displacement x ˜ that may be different from x. By inspection of a formal eigenfunction expansion of u, it is seen that

x x ˜ , (7)

if dominantly one eigenmode is excited in the spatial vibration pattern. The dominant biological mode is the fundamental mode, which is most often the excited vibrational mode of elastic structures relevant to auditory processes. The second assumption is that x x ˜ can be utilized as x = x ˜ . That is, only one mode of the elastic structure is excited dominantly.

Effectively, the partial differential equation that would be needed for a full description of one constituent of the tympanic structure can be written in the form of a harmonic oscillator equation,

m x ¨ + k x = F full . (8)

To derive the two-spring model we will use the following argument: Since the spherical shell is supported elastically by a ring underneath, we only need to store local vibrations of the shell surface in x s .

The displacement of the shell as a whole along the symmetry axis of the ring is equal to the displacement x r of the ring. With phenomenological damping coefficients γ s and γ r average displacements x s and x r in the sense of the piston approximation satisfy damped harmonic oscillator equations with two driving source terms F s , full and F r , full

m s x ¨ s + 2 m s γ s x ˙ s + k s x s = F s , full ( t ; x s , x r ) , (9)

m r x ¨ r + 2 m r γ r x ˙ r + k r x r = F r , full ( t ; x s , x r ) . (10)

We now need to specify the driving forces F s , full and F r , full in more detail. First, we assume that an external pressure signal p hits the spherical shell homogeneously. Put simple, p transfers a uniform momentum to the shell’s surface and thus makes the shell move up and down without deformation of its equilibrium shape. By definition of x s , x r , this corresponds to an external force applied to x r because x s shall only store local shell vibrations but no global piston modes. Since the ring itself is elastic, it will vibrate and thus lead to a small bending of the shell [15] .

For the elastic interaction between the two constituent elastic structures, we use the Hooke approximation. We regard the constituents as being coupled by a massless, elastic spring with spring constant k 0 : F s , full = k 0 x r and F r , full = F + k 0 x s where k 0 is treated as a fitting parameter. A rough estimate concerning the ratio k 0 / k r can be found in the Appendix. The other parameters m s , m r , γ s , γ r , k s , k r can be determined from experimental measurements. The above choice of source terms defines the two-spring model. Its underlying physical picture is summarized for the reader’s convenience in Figure 2, on the right.

2.2. Two-Spring Model

The deliberations of the previous paragraph are summarized by the governing equations of our two-spring model,

m s x ¨ s + 2 m s γ s x ˙ s + k s x s = k 0 x r , (11)

m r x ¨ r + 2 m r γ r x ˙ s + k r x r = k 0 x s + F . (12)

We assume k 0 < min { k s , k 0 } . This means that the coupling through the massless spring does not dominate the spring constants corresponding to eigenmodes of the individual constituents.

Since the tympanic structure is supposed to be at rest in the absence of an incident acoustic pressure signal on the shell, the initial conditions at t = 0 , i.e., when a signal is about to hit the shell, are fully homogeneous:

x s ( t = 0 ) = 0 = x r ( t = 0 ) and x ˙ r ( t = 0 ) = 0 = x ˙ s ( t = 0 ) .

Let us denote the surface of the shell by Σ ( t ) . An incident pressure wave is described as p ( t , x ) where x are the coordinates in the surrounding space. It becomes a vector quantity when we multiply it by the direction of wave propagation, viz., the normal to the wave front. Then the scalar net force exerted on the shell is given by

F ( t ) = ( Σ d S ( Σ ) , p 3 ) ( t ) . (13)

In practice, a time harmonic signal F ( t ) p 0 Σ 0 exp ( i ω t + i ϕ ) proves to give decent results for an incident low-frequency plane-wave-like signal with uniform scalar amplitude p 0 due to a sourced at a position in the acoustic far field. Σ 0 is the area of the equilibrium shape of the shell constituent, ω denotes the frequency carried by the acoustic signal, and ϕ symbolizes where to place an additive phase to account for directional information stored in the signal.

2.3. Solution of the Two-Spring Model

We solve the two-spring model by reduction to a first-order system. Let us define y s = x ˙ s and y r = x ˙ r . With a r F / m r , the first-order system reads

d d t ( x s y s x r y r ) = ( 0 0 0 a r ) + M ( x s y s x r y r ) where M ( 0 1 0 0 ω s 2 2 γ s ω d , s 2 0 0 0 0 1 ω d , r 2 0 ω r 2 2 γ r ) . (14)

We have employed the following definitions of squares of reduced frequencies: ω r 2 = k r / m r , ω s 2 = k s / m s and ω d , s 2 = k 0 / m s as well as ω d , r = k 0 / m r . The system is linear so that we can invoke the variation-of-constants or Duhamel formula [18] . We denote by [ A ] k l the ( k , l ) -entry of a matrix A . Conveniently, the initial conditions to be satisfied by x r , x s are homogeneous so that the variation of constants method yields

x s ( t ) = 0 t d τ [ exp ( ( t τ ) M ) ] 14 a r ( τ ) , (15)

x r ( t ) = 0 t d τ [ exp ( ( t τ ) M ) ] 34 a r ( τ ) (16)

2.4. Solution to the Eigenvalue Problem

The matrix exponential requires the solution of the eigenvalue problem for the matrix M and, in particular, knowledge of the zeros of the characteristic polynomial χ M ( ω ) det ( i ω 1 4 + M ) . More precisely, we need the roots of the following quartic equation,

0 = ( ω r 2 ω s 2 ω d , r 2 ω d , s 2 ) + 2 i ω ( ω r 2 γ s + γ r ω s 2 ) + ω 2 ( 4 γ r γ s ω r 2 ω s 2 ) 2 i ω 3 ( γ r + γ s ) + ω 4 d + c ω + b ω 2 + a ω 3 + ω 4 . (17)

The above equation can be analytically solved for ω by combining the methods of Cardano and Ferrari [19] [20] for the solution of a generic cubic (Cardano) and quartic (Ferrari) polynomial equations by algebraic completeness of . By a somewhat lengthy but straightforward calculations we can specify the four, possibly multiple, roots ω k ’s for k { 1,2,3,4 } in terms of

ω s 2 , ω d , s 2 , ω r 2 , ω d , r 2 , γ s , γ r :

ω 1 = [ i ( γ s + γ r ) p + 2 Ω 2 + p + 2 Ω 2 + 4 ( q 2 p + 2 Ω 2 + p + Ω 2 ) ] / 2 , (18)

ω 2 = [ i ( γ s + γ r ) p + 2 Ω 2 p + 2 Ω 2 + 4 ( q 2 p + 2 Ω 2 + p + Ω 2 ) ] / 2 , (19)

ω 3 = [ i ( γ s + γ r ) + p + 2 Ω 2 + p + 2 Ω 2 + 4 ( q 2 p + 2 Ω 2 + p + Ω 2 ) ] / 2 , (20)

ω 4 = [ i ( γ s + γ r ) + p + 2 Ω 2 p + 2 Ω 2 + 4 ( q 2 p + 2 Ω 2 + p + Ω 2 ) ] / 2 (21)

Ω 2 = 5 p 6 + 27 ( 125 p 3 108 5 p ( p 2 r ) 3 + 4 p 2 4 p r q 2 8 ) + 729 ( 125 p 3 108 5 p ( p 2 r ) 3 + 4 p 2 4 p r q 2 8 ) 2 + 108 ( 25 p 2 12 + 2 ( p 2 r ) ) 54 3 + 27 ( 125 p 3 108 5 p ( p 2 r ) 3 + 4 p 2 4 p r q 2 8 ) 729 ( 125 p 3 108 5 p ( p 2 r ) 3 + 4 p 2 4 p r q 2 8 ) 2 + 108 ( 25 p 2 12 + 2 ( p 2 r ) ) 54 3 (22)

p = 3 i ( γ r + γ s ) 3 4 γ r γ s ω r 2 ω s 2 (23)

q = i ( ( γ r + γ s ) ( 4 γ r γ s + ω r 2 + ω s 2 ) 2 ( ω r 2 γ s + γ r ω s 2 ) + ( γ r + γ s ) 3 ) , (24)

r = 3 ( γ s + γ r ) 4 16 + ( γ s + γ r ) 2 ( 4 γ r γ s + ω s 2 + ω r 2 ) 4 + ( γ s + γ r ) ( γ r ω s 2 + γ s ω r 2 ) + ( ω s 2 ω r 2 ω d , s 2 ω d , r 2 ) . (25)

2.5. Damping

It is interesting to note that damping enters directly through the average ( γ s + γ r ) / 2 of the individual damping coefficients γ r , γ s as specified in Equation (11) and Equation (12). Under physical assumptions regarding the measurement parameters, k s , k r , γ s , γ r , m s , m r and the fit parameter k 0 , we may assume that all four roots (18), (19), (20), and (21) are pairwise distinct. By a small variation of the parameters within the range of measurement uncertainty, this assumption can be fulfilled trivially. We denote by Z the resulting set of the four solutions ω to the eigenvalue problem of the matrix M .

2.6. Calculation of the Matrix Exponentials

To obtain x s and x r from (15) and (16), we need to calculate the matrix exponential. The calculation is performed with the aid of Mathematica.

We find

x s ( t ) = ω Z 0 t d τ G ω ( t τ ) a r ( τ ) , (26)

x r ( t ) = ω Z 0 t d τ H ω ( t τ ) a r ( τ ) (27)

where the sum runs over the four solutions ω 1 , ω 2 , ω 3 , ω 4 obtained in the previous paragraph. The kernels G ω ( t τ ) and H ω ( t τ ) are obtained from the eigenmode solutions upon diagonalization of M and defined such that

[ e ( t τ ) M ] 14 = ω Z G ω ( t τ ) , (28)

[ e ( t τ ) M ] 34 = ω Z H ω ( t τ ) . (29)

Assuming regularity of the input force F and thus of a r and by the assumption of distinctness of the ω k Z , the sum and the integral in the solution formulas (26) and (27) have been interchanged. Using the computational methods summarized above, we find the convolution kernels G ω and H ω to be given by

G ω ( t τ ) = ω d , r 2 2 e i ω ( t τ ) 2 i ω 3 3 ω 2 γ r 3 ω 2 γ s + 4 i ω γ r γ s + i ω ω r 2 + i ω ω s 2 + ω r 2 γ s + γ r ω s 2 (30)

H ω ( t τ ) = 1 2 ω 2 e i ω t τ + 2 i ω γ r e i ω t τ + ω 2 e i ω t τ 2 i ω 3 3 ω 2 γ r 3 ω 2 γ s + 4 i ω γ r γ s + i ω ω r 2 + i ω ω s 2 + ω r 2 γ s + γ r ω s 2 . (31)

Needless to say that although the model is simple and captures the physics of a combined movement of the spherical shell membrane due to bending and the elastic deformation of the supporting ring, the solution is somewhat cumbersome.

3. Numerics & Physical Discussion

3.1. Numerical Results

For the numerical simulation we use a time scale with units ω r 1 and a length scale of units F 0 / ( m r ω r 2 ) . The remaining parameters are fixed by γ r = 0.5 ω r , γ s = 0.7 ω r , ω s 2 = 0.8 ω r 2 and m s = m r . The fit constant k 0 is chosen such that ω d , s 2 = ω d , r 2 = 0.25 ω r 2 and ω = 0.9 ω r . In the top row of Figure 3, one sees that the elastic ring vibrates at a larger amplitude than the amplitude of bending of the spherical shell for both the real part (left) and the imaginary part (right) of the solution for the input exp ( i 0.8 t ) in the unit system specified above.

The color code employed in Figure 3 is as follows for the ring’s vertical elastic displacement x r . The real part [ x r ] in the top left plot and the imaginary part [ x r ] in the top right plot are shown as solid blue curves. As for the shell’s bending displacement x s , the following color code is employed in Figure 3. The real part [ x s ] in the top left plot and the imaginary part [ x s ] in the top right plot are shown is solid orange curves. Varying the parameter k 0 in the model simulation, it is found that the smaller k 0 , the smaller the amplitude of x s .

For the same choice of parameters, we have plotted the ratio

| [ x s ] | / ( | [ x s ] | + | [ x r ] | )

in the bottom-left plot of Figure 3 and

| [ x s ] | / ( | [ x s ] | + | [ x r ] | )

in the bottom-right plot of Figure 3. By definition, both ratios are smaller than unity. This is indicated by the orange lines in the two bottom plots. The ratios

Figure 3. Top-left plot: [ x r ] (blue) and [ x s ] (orange) in units of F 0 / ( m r ω r 2 ) over time t [ 0,15 ] in units of ω r 1 . Bottom-left: Absolute displacement over sum of absolute displacements | [ x s ] | / ( | [ x s ] | + | [ x r ] | ) as a function of t [ 0,15 ] measured in units of ω r 1 Top-and bottom-right plots: Analogous plots as in the left column with the corresponding imaginary parts used instead of the real parts.

quantify how important the bending of the plate is in comparison to the vertical oscillation of the ring at a given time t. The peaks in the bottom row plots in Figure 3 result from distinct zeros of [ x s ] , [ x r ] and [ x s ] , [ x r ] . Yet, and most of the time, the bending of the spherical shell is small compared to the vertical displacement of the ring.

The amplitude ratios max [ x s ] / max [ x r ] and max [ x s ] / max [ x r ] have been found by looking for the maximum amplitudes in t [ 4,15 ] after the periodic behavior stabilizes in the range t 4 because of relaxation of the system (by damping) into the quasi-stationary state. We found

max [ x s ] / max [ x r ] 0.2 max [ x s ] / max [ x r ] ,

rounded mathematically to one counting digit. This agrees with the reasoning of the introduction that the elastic ring is the preferred constituent of the tympanic structure of multitudinous archosaurs to generate the fundamental frequency that is processed further by ICE. As for the extracolumella that is wired to the cochlea, the relative smallness of the bending of the shell compared to the elastic stretching of the ring implies that the its deflection from equilibrium, x s + x r , can be approximated as x r for sufficiently small k 0 compared to k r and k s .

3.2. Physical Discussion

The quantity of keen interest to model builders in auditory research is the velocity of the tympanic membrane, or, more generally, of the tympanic structure [13] . While for geckos and frogs the interaural cavity can be and has been approximated by a cylinderical acoustic waveguide of equivalent volume [3] [4] [6] , the archosaur counterpart has the topology of a, sometimes even higher-genus, torus [21] . In view of ongoing analytic approaches to the study of hearing in “icy” animals, the geometric complexity as well as the number of equations required for mathematical description both necessitate approximations for efficient solutions. The piston approximation in conjunction with approximate modeling of the entire interaural cavity as effectively one-dimensional wave guide form the base of the preferential toolkit.

For comparison with biological data, quantities of the structure

ILD = 20 log 10 | u ˙ L u ˙ R | (32)

have proved useful [3] [4] . Such quantities are inspired by the definition of the sound pressure level in acoustics. Via no-penetration boundary conditions ρ air t ( t u ) = ( n ^ p ) at the eardrums, these quantities compare the acoustic pressure at the left to the one at the right ear. In Equation (32), u ˙ L and u ˙ R denote the (total) speed of the left and right tympanic membrane. In the new case of a tympanic structure with more than one constituent, a similar measure can be obtained by considering the total displacement of the tympanic structure, denoted by U L and U R . In the case of the present two-spring model, we identify U ˙ L = x ˙ s , L + x ˙ r , L and U ˙ R = x ˙ s , R + x ˙ r , R for the total speed at the left and right tympanic structure. The substitutions u L U L and u R U R in (32) result in

ILD = 20 log 10 | x ˙ s , L + x ˙ r , L x ˙ s , R + x ˙ r , R | . (33)

The primary objective of this work is to provide auditory research with a practically applicable model for the tympanic structure of numerous archosaurs. We discuss a simpler treatment of the two-spring model in Section 4.1.

4. Approximate Solution

4.1. Iterative Solution of the Two-Spring Model

Since the constituents of the tympanic structure―see (11) and (12)―are damped, an iterative solution is to be favored instead of the lengthy closed form solution derived in the previous paragraph. By iterative solution, we mean a decoupling of the equations of motion for the full tympanic structure into the individual equations of motion for its constituents using the assumption k 0 < min { k r , k s } .

Because of damping, the impact of higher-order coupled contributions to the displacement of either constituent from the other constituent will decrease in time. We refer to Hassani [18] for a mathematical introduction to operators as used in physics. For the sake of notational brevity, we introduce the following constant-coefficient differential operators

D s 2 m s d t 2 + 2 m s γ s d t + k s , (34)

D r 2 m r d t 2 + 2 m r γ r d t + k r . (35)

They are defined on the space H = C 2 ( 0 + ) . Moreover, we require the driving force F to be smooth. Both operators are accompanied by homogeneous initial conditions. We have [ D r 2 , D s 2 ] = 0 on H C 4 ( 0 + ) . That is, the operators defined in (34) and (35) commute for high enough regular functions. The governing equations of the two-spring model read in operator notation

D s 2 [ x s ] = k 0 x r & D r 2 [ x r ] = k 0 x s + F . (36)

Using one of Schwarz’ representation theorems, we can represent the inverse operator by a convolution-type integral [18] . Let C 0 ( 0 + ) and t 0 . We then find

D s 2 [ ] ( t ) = 0 t d τ e γ s ( t τ ) sin ( ω s 2 γ s 2 ( t τ ) ) m s ω s 2 γ s 2 ( τ ) , (37)

D r 2 [ ] ( t ) = 0 t d τ e γ r ( t τ ) sin ( ω r 2 γ r 2 ( t τ ) ) m r ω r 2 γ r 2 ( τ ) . (38)

By direct computation, it is seen that D r 2 D r 2 = 1 = D r 2 D r 2 on H and analogously for D s 2 and D s 2 .

Next, the smoothness assumption F C ( 0 + ) allows us to confine the search to smooth solutions x s , x r H C ( 0 + ) . If we are able to find such x s and x r , the existence and uniqueness theorems from the theory of ordinary differential equations guarantee that we have also found the only solution to the inhomogeneous initial value problem. Using the inverse operators (37) and (38), we recast (36) into the equivalent system of integral equations

x s = k 0 D s 2 [ x r ] & x r = k 0 D r 2 [ x s ] + D r 2 [ F ] (39)

The Banach fixed-point theorem entails a construction to obtain the fixed-point by iteration. In order to satisfy the model’s initial conditions, we choose ( x s , x r ) = ( 0,0 ) H . Let t > 0 and define for all n by

( t 0 = 0 , t 1 , , t n = t ) ( 0 + ) n

with t i < t i + 1 for 0 < i < n 1 a subdivision of [ 0, t ] . We introduce the notation for subdivisions ( t 0 , t 1 , t 2 ) ( 0 + ) 3 ,

( t 0 , t 1 , t 2 , , t 2 k 1 , t 2 k ) ( 0 + ) 2 k + 1 , k > 1

with the properties declared above

A t 1 , t 2 [ ( t 0 ) ] ( t 2 ) D s 2 [ D r 2 [ ( t 0 ) ] ( t 1 ) ] ( t 2 ) , A t 1 , , t 2 k A t 1 , , t 2 k 2 A t 1 , t 2 . (40)

The indices have been inserted in the definition of the A operators to make the nesting of integrals explicit. They indicate the integration variables in the nested integrals.

For t > 0 , we find the following series representation for the components of ( x s , x r )

x s ( t ) = ( m = 1 k 0 m A t 1 , , t 2 m [ F ( t 0 ) ] ) ( t ) , (41)

x r ( t ) = ( m = 0 k 0 m D r 2 [ A t 1 , , t 2 m [ F ( t 0 ) ] ] ) ( t ) . (42)

The above expression is not a Neumann series because of the appearance of the nested integrals. It cannot be summed up to a closed-form solution. Recalling the fixed-point argument from above, we see that the maximum norm of the series converges. The regularity properties of x s , x r are the result of smoothness of the input F and the fact that the operators (37) and (38) have smooth integration kernels and thus do not decrease the regularity of their arguments. Equations (41) and (42) can be applied in practice by truncating the series. In the next paragraph, we compare truncations of the series solutions (41) and (42) with the exact solutions of the two-spring model for a given parameter set.

4.2. Discussion of the Applicability in Auditory Research

We state the first two iterations. In the notation of the previous paragraph, this corresponds to x s ( N ) , x r ( N ) with N { 1,2 } ,

x r ( 1 ) ( t ) = 0 t d τ e γ r ( t τ ) sin ( ω r 2 γ r 2 ( t τ ) ) m r ω r 2 γ r 2 F ( τ ) , (43)

x s ( 1 ) ( t ) = ω d , r 2 0 t d τ e γ s ( t τ ) sin ( ω s 2 γ s 2 ( t τ ) ) m s ω s 2 γ s 2 [ 0 τ d τ e γ r ( τ τ ) sin ( ω r 2 γ r 2 ( τ τ ) ) ω r 2 γ r 2 F ( τ ) ] , (44)

x r ( 2 ) ( t ) = 0 t d τ e γ r ( t τ ) sin ( ω r 2 γ r 2 ( t τ ) ) m r ω r 2 γ r 2 F ( τ ) + ω d , r 2 ω d , s 2 0 t d τ e γ r ( t τ ) sin ( ω r 2 γ r 2 ( t τ ) ) m r ω r 2 γ r 2 { 0 τ d τ e γ s ( τ τ ) sin ( ω s 2 γ s 2 ( τ τ ) ) ω s 2 γ s 2 [ 0 τ d τ e γ r ( τ τ ) sin ( ω r 2 γ r 2 ( τ τ ) ) ω r 2 γ r 2 F ( τ ) ] } , (45)

x s ( 2 ) ( t ) = ω d , r 2 0 t d τ e γ s ( t τ ) sin ( ω s 2 γ s 2 ( t τ ) ) m s ω s 2 γ s 2 [ 0 τ d τ e γ r ( τ τ ) sin ( ω r 2 γ r 2 ( τ τ ) ) ω r 2 γ r 2 F ( τ ) ] + ω d , s 2 ω d , r 4 0 t d τ e γ s ( t τ ) sin ( ω s 2 γ s 2 ( t τ ) ) m s ω s 2 γ s 2 { 0 τ d τ e γ r ( τ τ ) sin ( ω r 2 γ r 2 ( τ τ ) ) ω r 2 γ r 2 [ 0 τ d τ e γ s ( τ τ ) sin ( ω s 2 γ s 2 ( τ τ ) ) ω s 2 γ s 2

( 0 τ d τ e γ r ( τ τ ) sin ( ω r 2 γ r 2 ( τ τ ) ) ω r 2 γ r 2 F ( τ ) ) ] } . (46)

From the phenomenological point of view, the “strength” of coupling can be determined from the ratios ω d , s 2 / ω s 2 = k 0 / k s or ω d , r 2 / ω r 2 = k 0 / k r for 2-constituent tympanic structures. By the assumptions regarding k 0 , both ratios are smaller than unity. A suitable smallness parameter is defined through ϵ = max { k 0 / k s , k 0 / k r } . By mathematical induction it follows

| x s ( N ) x s ( N 1 ) | / | x s ( N 1 ) | | x r ( N ) x r ( N 1 ) | / | x s ( N 1 ) | ϵ 2 .

The series solution gradually builds up the coupling between the two constituents of the tympanic structure.

4.3. Numerical Results

In order to assess the accuracy of the iterative solution developed in this section, we simulate the Equation (43) and Equation (44), respectively (45) and (46), for both the parameter choice and the input model used for the simulation of the full two-spring model presented the previous section. The numerical results for (43) and (44) are shown in Figure 4. The axis units therein have been chosen in accordance with Figure 3 and the unit system introduced at the beginning of the previous section. The numerical results for (45) and (46) are shown in Figure 5. In each of Figure 4 and Figure 5, the left column depicts results for the real parts [ x r ( 1 ) ] and [ x s ( 1 ) ] , respectively [ x r ( 2 ) ] and [ x s ( 2 ) ] .

The right columns of Figure 4 and Figure 5 show the results for [ x r ( 1 ) ] and [ x s ( 1 ) ] , respectively [ x r ( 2 ) ] and [ x s ( 2 ) ] . In the top row of Figure 4, [ x r ( 1 ) ] , I m [ x r ( 1 ) ] is shown as solid blue curves, while [ x s ( 1 ) ] , I m [ x s ( 1 ) ] is represented by the solid orange curves. The same color code has been employed in the top row of Figure 5 with x r ( 2 ) and x s ( 2 ) used in place of x r ( 1 ) and x s ( 1 ) .

The bottom rows of the two figures are devoted to quantifying the absolute deviations from the full solution obtained by means of the simulation presented in the previous section. In the two bottom plots of Figure 4, the differences | [ x ] r ( 1 ) | | [ x ] r | and | [ x ] s ( 1 ) | | [ x ] s | are shown as the solid blue and solid orange curve on the left and the differences | [ x ] r ( 1 ) | | [ x ] r | and | [ x ] s ( 1 ) | | [ x ] s | are shown as the solid blue and solid orange curve on the right. Similarly, the differences | [ x ] r ( 1 ) | | [ x ] r | and | [ x ] s ( 1 ) | | [ x ] s | are shown as the solid blue and solid orange curve in the left plot and the differences | [ x ] r ( 1 ) | | [ x ] r | and | [ x ] s ( 1 ) | | [ x ] s | are shown as the solid blue and solid orange curve in the right plot of the bottom row in Figure 5.

At a first sight, the top rows of Figure 4 and Figure 5 look like their “exact” pendants in the top row of Figure 5. We see, however, that the deviations of the moduli of the real and imaginary part of the first iterates, i.e., (43) and (44), from the moduli of the real and imaginary part of the exact solution of the two-spring model do not equilibrate fast into the quasi-stationary state behavior; cf. the bottom row of Figure 4.

In contrast, the bottom row of Figure 5 shows that, for the specified parameter values, the moduli of the real and imaginary part of the exact solution of the two-spring model and the absolute values (45) and (46) of the real and imaginary part of the second iterates match well, as compared with the amplitude of the exact solution, after equilibration into the quasi-stationary state. For t 6 ω r 1 , the

Figure 4. Top-left: Approximate solution [ x r ( 1 ) ] (blue) and [ x s ( 1 ) ] (orange) in units of F 0 / ( m r ω r 2 ) after the first iteration as a function of t [ 0,15 ] in units of ω r 1 . Bottom-left: Comparison of the solution for the full two-spring model with the approximate solutions: | [ x r ( 1 ) ] | | [ x r ] | (blue) and | [ x s ( 1 ) ] | | [ x s ] | (orange) in units of F 0 / ( m r ω r 2 ) as a function of t [ 0,15 ] in units of ω r 1 Top- and bottom right: Analogous plots as in the left column with the corresponding imaginary parts used in place of the real parts.

Figure 5. Top-left: Approximate solution [ x r ( 2 ) ] (blue) and [ x s ( 2 ) ] (orange) in units of F 0 / ( m r ω r 2 ) after the second iteration as a function of t [ 0,15 ] in units of ω r 1 . Bottom-left: Comparison of the solution for the full two-spring model with the approximate solutions: | [ x r ( 2 ) ] | | [ x r ] | (blue) and | [ x s ( 2 ) ] | | [ x s ] | (orange) in units of F 0 / ( m r ω r 2 ) as a function of t [ 0,15 ] in units of ω r 1 . Top- and bottom-right: Analogous plots as in the left row with the corresponding imaginary parts used in place of the real parts.

relative deviation compared to the amplitudes is below 1% for both real and imaginary part for both the ring and the shell displacement. This is lower than the relative deviation compared to the amplitude of the oscillation of the real and imaginary part of x r and x s found for the first iterates from the bottom row of Figure 4. In this case, we find the ratio in question being under 6% for t 6 ω r 1 .

5. Summary

For the first time in auditory research, we have articulated a physically implementable definition of tympanic structures that generalizes the notion of a tympanic membrane. In the terminology introduced in the main body of the present article, a tympanic membrane is seen to be a special case of the wider notion of “tympanic structure”. The case of archosaurs as “icy” animals exemplified the necessity to define a more general notion in order to model the tympanic responses in these animals in a physically meaningful way. Favoring applicability of the ICE model in auditory research over mathematical complexity, the two-spring model has been derived from the piston approximation applied to the elastic constituents of the 2-constituent tympanic structure as observed in numerous archosaurs. Two damped harmonic oscillators, interacting via a Hooke-elastic coupling, were found to account for the combined motion of the tympanic structure. That is, the displacement of the spherical shell membrane of the animals without local bending and bending of the shell.

By numerical simulations with parameters chosen from the expected parameter range, the biologically motivated hypothesis that the overall displacement of the extracolumella is mainly due to uniform, piston-like displacement of the shell as a whole without large contributions due to local bending, could be supported. Only a maximum of 20% of the extracolumella displacement’s amplitude was found to be on the part of shell bending in the two-spring model with the chosen simulation parameters.

Because of the length of the analytic solution formulas, a simpler iterative treatment has been developed. Treating the two-spring model as a set of two damped, suitably weakly coupled harmonic oscillators, a Banach fixed-point iteration scheme has been derived. Simulations of the first two non-trivial iterates for the model solutions were performed using the same parameters as in the exact solution of the two-spring model. For sufficiently weak coupling of the two constituents of the tympanic structure, the second iterates were found to be sufficiently accurate and efficient-to-calculate at the same time so as to be favored over the exact solution, if the computational power available is limited. The strength of the coupling is stored in the parameter ϵ = max { k 0 / k r , k 0 / k r } , which requires specification through up-to-now unavailable experimental data.

Acknowledgements

The authors thank Professor Daniel Rixen (Mechanical Engineering, TU Munich) for helpful suggestions regarding the literature on plates and shells, and Professor Bruce Young (A.T. Still University, Kirksville, MO) for insightful remarks concerning the anatomy of crocodilian eardrums. They apologize to the latter for biologically overly simplifying a mathematically highly complicated problem.

Appendix

The two-spring model is an effective model, a strong simplification of both biological and mathematical reality. In order to ensure the overall applicability of the model as well as the readability for a less mathematically inclined audience, a more complete, yet mathematically more advanced model, has been postponed to this appendix. Here spheroidal wave functions will appear at the end. Unfortunately, the state of the art in the numerics of spheroidal wave functions obstructs a detailed numerical evaluation while further analytic investigations are ruled out by the non-existence of closed form representations of spheroidal wave functions.

Geometry of the configuration

As explained in the main text, the two constituents of the tympanic structure in multitudinous archosaurs are (1) an elastic sectorial spherical shell and (2) an elastic circular symmetric ring. In spherical coordinates, the shell has radius R 0 , azimuthal opening angle θ 0 [ 0, π / 2 ] , and a polar opening angle ϕ [ β ,2 π β ] where 0 < β π . The breaking of the full rotational symmetry of the tympanic structure results from the presence of the extracolumella; cf. Figure 1. It is a bony structure placed in the sector ϕ [ 0 , β ) ( 2 π β , 2 π ] of the shell. The extracolumella effectively prevents vibrations on the shell surface itself. Yet, in the present model the shell is allowed to move up and down, the entire shell moving uniformly as a piston.

For the ensuing treatment, the extracolumella is assumed to be only connected to the shell but not to the underneath ring. The ring itself has radius R 0 and height H, extending in negative z-direction. We denote by Ω s 0 , Ω r 0 3 the equilibrium surfaces of shell and ring,

Ω s 0 = S R 0 2 ( 0 3 ) { ( x , y , z ) 3 | z 0 } , (47)

Ω r 0 = S R 0 1 ( 0 2 ) × [ 0 , H ] . (48)

For the elasto-dynamic analysis, it is useful to account for the small but finite thickness of the elastic constituents as well. Let us denote the corresponding compacta by Ω s d s and Ω r d r where d s , d r min { R 0 , H } stand for the thickness of the shell and the ring, respectively. We assume d s = d r and use d s as a symbol for the thickness of shell and ring. Explicitly, Ω r d s and Ω s d s are defined as

Ω s d s = Cl ( B R 0 + d s / 2 3 ( 0 3 ) \ B R 0 d s / 2 3 ( 0 3 ) ) { ( x , y , z ) 3 | z 0 } , (49)

Ω r d s = Cl ( B R 0 + d s / 2 2 ( 0 2 ) \ B R 0 d s / 2 2 ( 0 2 ) ) × [ 0 , H 0 ] (50)

where Cl ( M ) denotes the topological closure operation on the set M in brackets. Letting d s 0 tend to zero from above, we reproduce Ω s 0 and Ω r 0 as the maximal set contained in all the Ω s d s ’s and Ω r d s ’s, respectively. The interface between the two constituents is given the intersection Ω s d s Ω r d s Γ d s . It is localized in the z = 0 -plane, Γ d s 2 × { 0 } and is isometrically homeomorphic to the closed annulus A ¯ R 0 d s / 2 , R 0 + d s / 2 . If d s 0 tends to zero from above, the intersection of all the Γ s d s ’s gives Γ 0 S R 0 1 ( 0 2 ) × { 0 } .

Physical description of the problem

Each of the two constituents of the tympanic structure is an elastic material that satisfies approximately material isotropy and homogeneity assumptions. That is, shell and ring have respective Lamé constants ( λ s , μ s ) and ( λ r , μ r ) . Since, even in the absence of a physical cause for deformations, the shell plus extracolumella has the shape of a hemisphere, it is treated as self-supporting: If stimulated by an pressure signal incident from the z-direction, the shell undergoes elastic deformations only along the z-direction.

Navier-Cauchy equation

The continuum momentum balance equations simplify to linear partial differential equations if only small elastic deformations are to be modeled. For the shell and the ring, they read quite generally

ρ s t 2 u s = σ s + f s , (51)

ρ r t 2 u r = σ r + f r (52)

where ρ s , ρ r denote the constant mass per volume density of shell and ring material, respectively, f s and f r stand for an external force density applied to the shell and ring, and σ r and σ r denote the (symmetric) stress tensors. The Navier-Cauchy equations [23] are obtained as a consequence of the assumption of a linear stress-strain relation. The linearized strain tensors ϵ r = 1 / 2 ( ( u r ) + ( u r ) T ) and ϵ s = 1 / 2 ( ( u s ) + ( u s ) T ) for the ring and shell occur in the linear stress-strain relation as

( σ s ) i j = ( Λ ( s ) : : ϵ s ) i j and ( σ r ) i j = ( Λ ( r ) : : ϵ r ) i j .

The : : denotes double contraction in engineers’ notation. The material homogeneity and isotropy assumptions reduce the number of components from a priori 3 4 = 81 to 2 for both of the rank 4 tensors Λ ( r ) and Λ ( s ) . The stress tensors σ r and σ s are themselves symmetric. They can only depend linearly on the respective strain tensor and the trace of the respective strain tensors. A textbook calculation [23] demonstrates that the Equation (51) and Equation (52) reduce to the elasto-dynamic equations, also known as the Navier-Cauchy equations,

ρ s t 2 u s = ( λ s + 2 μ s ) ( u s ) μ s × × u s + f s , (53)

ρ r t 2 u r = ( λ r + 2 μ r ) ( u r ) μ r × × u r + f r . (54)

Relations of the Lamé constants to other elastic parameters

Let λ and μ indicate Lamé’s constants. The literature denotes μ occasionally by G and calls G the shear modulus. The expression λ + μ = G / ( 1 2 ν p ) defines the Poisson number ν and is smaller than 0.5. The Young’s modulus E satisfies E = 2 G ( 1 + ν ) = 4 μ 2 / ( 2 μ λ ) . The bulk modulus K relates to the Lamé constants via K = E / ( 3 ( 1 2 ν p ) ) = ( λ + μ ) ( 2 μ λ ) / ( 3 μ ) .

Model specification

The shell’s force-density drive is assumed to take place only in the z-direction due to the self-supporting property of the shell. The prototypical external stimulus used in auditory research is a pressure signal emanating from a point far away from the shell surface so that a plane is the more appropriate shape for the wave front than a shell. Since the propagation speed of sound times a few milliseconds still exceeds the geometric dimensions of the tympanic structure, arrival time differences on the shell surfaces are negligible. Thus the net force-per-volume density on the shell is spatially uniform and purely time-dependent, and only has a non-trivial z-component: f s = f z ( s ) ( t ) e ^ z .

Its relation to the acoustic pressure signal is established by considering the z-component of force-per-volume density to be the pressure difference along the thickness of the shell: f z ( s ) ( t ) = p ( t ) / d s . Indeed, integrating over the homogeneous membrane’s thickness d s , we reproduce the pressure difference p ( t ) as the force-per-surface density which is responsible for driving the system.

The ICE model has proven quite successful in auditory research [1] [12] . Due to the superposition principle, the investigation is confined to a pure tone, p ( t ) = p 0 exp ( i ω t ) with a fixed positive stimulus angular frequency ω > 0 . The ring only interacts with the shell but receives no external stimulation, i.e., f r = 0 . The ring-shell interaction via boundary conditions will be discussed later on.

Our focus is on the quasi-stationary state. That is, the tympanic structure is assessed at a time where tympana’s undulations generated previously have decayed as a result of damping. Mathematically, this amounts to working in frequency-domain space where linearity ensures that the tympana undulate with the same angular frequency ω as carried by the stimulus but in general with a phase shift for stimulus directions ¹ 0. The self-supporting assumption reduces the components of the displacement vector field u s from 3 to only 1, namely, to the component aligned in parallel to the external stimulus: u s = u s ( z ) e ^ z .

As mentioned before, the shell rests on top of the ring. If piston-like motion is excited, the entire ring vibrates as if it were a point mass under the influence of Newton’s law in vertical direction. Neglecting back-couplings of the ring’s displacement to spatially dependent vibrations on the shell, the ring’s displacement is aligned in parallel to u s so that u r = u z ( r ) ( t ) e ^ z . The configuration hints at using cylindrical coordinates for the spatial arguments of the functions u z ( s ) and u z ( r ) . Upon insertion of the vector operators in cylindrical coordinates [22] , (53) and (54) simplify,

ρ s t 2 u z ( s ) = ( λ s + μ s ) z 2 u z ( s ) + μ s Δ s u z ( s ) + p ( t ) / d s , (55)

ρ r t 2 u z ( r ) = ( λ r + μ r ) z 2 u z ( r ) + μ r Δ r u z ( r ) . (56)

Δ s and Δ r denote the Laplace-Beltrami operators in cylindrical coordinates on the three-dimensional equilibria volumes of shell and ring. That is, on Ω s d s and Ω r d s . The expression Δ v = ( ) × × can be utilized to convert (54) and (53) to hyperbolic differential equations akin to (56) and (55) featuring only the vector Laplacian and the composition of gradient and divergence as differential operators in spatial variables.

Compressional coupling

We now deal with appropriate elastic matching conditions at the interface of ring and shell, i.e., on Γ d s . Requiring validity of the differential Equation (56) and Equation (55) on Ω r d s Γ d s and Ω s d s Γ d s , respectively, ensures that u z ( s ) and u z ( r ) are C 2 on Γ d s as well. As a side comment this means that we would be allowed to extend the solutions to a larger but open domain, say, Ω r , + Ω r d s and Ω s , + Ω s d s , so that we can already assume the solutions to be sufficiently well-behaved.

We recall that the compression of an elastic medium generates an excess pressure inside that medium through p m = K m u m . The index m refers to the medium, whence m { r , s } here. Upon insertion of the divergence operator in cylindrical coordinates and due to the ansatz for u r and u s from the previous paragraph, the excess pressures p r and p s are given by

p r = K r z u z ( r ) (57)

p s = K s z u z ( s ) (58)

The regularity requirement discussed at the beginning of this paragraph permits the usage of these expressions on the interface Γ d s . The difference in material generates a pressure difference δ p = p r p s = K r z u z ( r ) + K s z u z ( s ) localized at z = 0 . One possible boundary condition would be to equate this compressional pressure difference and its partial derivative in z-direction to zero on Γ d s . In the configuration under consideration this is, however, not a good idea because the ring experiences yet another pressure on Γ s d s . Viz., the one generated by pressure signal p ( t ) , which pulls the shell as a mass point up and down. The detailed discussion of the consequences is the subject of the next paragraph.

Here we confine ourselves to introducing the compressional coupling conditions at the interface Γ s d s . In Newton’s second law, the gradient of the pressure difference δ p generates a force-per-volume. The boundary of shell and ring vibrate in such a way that the pressure gradients in positive and negative z-direction are compensated,

ρ s t 2 u z ( s ) = z δ p = K r z 2 u z ( r ) + K s z 2 u z ( s ) , (59)

ρ r t 2 u z ( r ) = z δ p = K r z 2 u z ( r ) K s z 2 u z ( s ) . (60)

The equations express mathematically that the two constituent structures can be compressed at the interface. The resulting difference in excess pressures, δ p , acts as force drive for the shell and the ring at the interface in the z = 0 plane. The bulk moduli K r and K s are material constants that can be determined from the Lamé constants ( λ r , μ r ) and ( λ s , μ s ) for the ring and shell, respectively, see the paragraph entitled “Relations of the Lamé constants to other elastic parameters”. In order to deal with the piston mode where the total of shell and extracolumella is driven as one point mass by the incident pressure signal needs to be considered in the above equations.

Force-collector model

Let us suppose that the shell is placed on an incompressible iron floor instead of an elastic ring. Since the shell is self-supporting and the pressure uniform, only forces in z-direction need to be considered. By Newton’s third law, the iron floor enacts a excess normal force on the annular interface region of floor and shell as a reaction to the total of excess pressure driving the membrane. The resulting reaction force equates to the reaction pressure p int times the area of the interface A int . Since the shell obviously does not oscillate up and down through the “iron” force, the reaction force needs to compensate the total force generated by the incident pressure p on the area of the shell plus extracolumella as seen from a point on the z-axis, A hit = π R 0 2 . We note A hit A hemisphere = 2 π R 0 2 because radial oscillations of the sphere are prohibited by the self-supporting assumption. The interface region Γ d s between shell and floor is annular and has the area

A int = π ( R 0 + d s / 2 ) 2 π ( R 0 d s / 2 ) 2 2 π R 0 d s

neglecting contributions of order ~ d s 2 / R 0 2 relative to A hit . The global force balance according to Newton’s third law (actio = -reactio) translates into

p int = p ext A hit / A int = p ext R 0 / ( 2 d s ) .

The difference between the iron floor and the elastic ring is that the latter experiences the interface pressure as an additional pressure. We observe that (59) needs no correction, if the piston mode is subtracted from u z ( s ) since u z ( s ) only quantifies dislocations generated by the compressional coupling between shell and ring.

One technical issue needs to be fixed. In the shell model (55) the pressure difference along the shell interface enters the force-per-volume balance rather than only an incident pressure gradient. The shell collects the pressure difference, the physical motivation for the nomenclature “force-collector model”. It needs to be shown that the interface region Γ d s is subject to the collected pressure difference. Namely, the force-collector model demonstrates that an additional force-per-volume density acts on the ring at the boundary namely

f ext = R 0 / ( 2 d s ) p ( t ) / d s .

Shell and ring extend in radial direction by d s . Let us replace the fields u z ( s ) and u z ( r ) in (59) and (60) by their respective average along a ray segment of length d s in the { z = 0 } plane. After multiplication by d s , Equation (59) and Equation (60) read

σ s t 2 u z , avg ( s ) ( t , z = 0 , ϕ ) = ( K s d s ) z 2 u z , avg ( s ) ( t , z = 0 , ϕ ) ( K r d s ) z 2 u z , avg ( r ) ( t , z = 0 , ϕ ) , (61)

σ r t 2 u z , avg ( r ) ( t , z = 0 , ϕ ) = ( K r d s ) z 2 u z , avg ( r ) ( t , z = 0 , ϕ ) ( K s d s ) z 2 u z , avg ( s ) ( t , z = 0 , ϕ ) . (62)

The quantities σ r and σ s are effective surface mass densities defined through σ r ρ r d s and σ s ρ s d s . Given that d s R 0 , H , the approximation is sensible. Its advantage lies in the fact that p int as specified at the beginning of the paragraph can be included in these boundary conditions in the spirit of the “force collector” argument,

σ s t 2 u z , avg ( s ) ( t , z = 0 , ϕ ) = ( K s d s ) z 2 u z , avg ( s ) ( t , z = 0 , ϕ ) ( K r d s ) z 2 u z , avg ( r ) ( t , z = 0 , ϕ ) , (63)

σ r t 2 u z , avg ( r ) ( t , z = 0 , ϕ ) = ( K r d s ) z 2 u z , avg ( r ) ( t , z = 0 , ϕ ) ( K s d s ) z 2 u z , avg ( s ) ( t , z = 0 , ϕ ) + a p ( t ) . (64)

The dimensionless quantity a R 0 / ( 2 d s ) is called amplification factor. Upon division by d s , the statement that also the pressure difference along the interface is collected holds true in an approximate sense. Namely, we need to replace the fields u z ( r ) and u z ( s ) by the averages over the thickness of the annulus.

Effectively, we are left with fields on the interface Γ 0 rather than Γ d s . Indeed, d s H , R 0 already hints at using Ω s 0 and Ω r 0 in place of Ω s d s and Ω r d s . We need to answer the question of whether or not to keep the source term in (55). In order to address this issue, we need to asses the ratio of maximum amplitudes of shell and ring. At the sphere’s top at ( z = H 0 ) , the shell experiences a pressure as σ s t 2 u z , avg ( s ) ( t , z = H 0 ) = p ( t ) . Neglecting spatial variations on the interface for the estimate, (64) reduces to

σ r t 2 u z , avg ( r ) ( t , z = 0 , ϕ ) = a p ( t ) .

Since p ( t ) p 0 exp ( i ω t ) , the quotient becomes

u z , avg ( s ) ( t , z = H 0 ) u z , avg ( r ) ( t , z = 0, ϕ ) = σ s σ r 1 a = ρ s ρ r 2 d s R 0 1 (65)

in the quasi-stationary state, i.e., only the solution to the inhomogeneous problem is considered. The “ ” awaits some justification.

Typically, ρ s O ( ρ r ) because the material of the shell and ring are biological tissue. That is, mostly water. However, the thickness of the membrane d s 10 5 m R 0 exceeds R 0 0.5 × 10 2 m by 2 to 4 orders of magnitude! Compared to the elastic deformation of the ring which experiences an amplification a p 0 of the pressure collected by the shell, the shell only a pressure of amplitude p 0 . The amplification factor a contains purely geometrical information, a = R 0 / ( 2 d s ) = A hit / A int . The estimate of the orders of magnitude shows that the maximum amplitude of the piston shell is small compared to the amplitude of the ring compression. The “up and down” oscillation of the shell plus extracolumella system is thus dominantly due to the compressional oscillation of the elastic ring as a response to the total pressure. It therefore makes sense to drop the source term in (55) and account for the external stimulus in (64) instead.

The limit of thin shells and rings, d s 0 +

Since we are only interested in displacement averaged over the shell’s thickness d s , the variables r and z are no longer independent but subject to the constraint r = r ( z ) = R 0 2 z 2 . The clamping due to the presence of the extracolumella is unproblematic and only settles ϕ [ β ,2 π β ] . The constraint r 2 + z 2 = R 0 2 with z 0 is a cue to try transforming the differential operator of the spatial variables in (55) into spherical coordinates. The problem, however, is that λ s + μ s 0 . The case λ s + μ s = 0 corresponds to a negative Poisson number, more precisely ν = 1 . The relation between the shell’s Lamé constants λ s and μ s on the one hand and the Young’s modulus E and the Poisson number ν on the other hand becomes ill-defined. Furthermore, the Poisson number is typically between 0 and 1/2. Biological tissue does not increase its volume in such a way that it gets thicker under application of an external pulling force. Rather, its volume stays constant or increases while the shell becomes thinner [23] .

In view of the above issues, a different method to localize the shell on r 2 + z 2 = R 0 2 for z 0 is what we look for. One possibility starts with redefining z z where > 0 is a suitable re-scaling factor involving the Lamé constants. It is then chosen so that the problematic partial differential operator with spatial variables in (55) equals the Laplacian in the new coordinates. is then determined to be

2 + λ s μ s > 1. (66)

In the new coordinates ( r , ϕ , z ) , the physical hemisphere defined through r 2 + z 2 = R 0 2 with z 0 is half of an oblate spheroidal hemisphere,

r 2 R 0 2 + ( z ) 2 R 0 2 2 = 1. (67)

Implementing the deliberations outlined above, (55) reads in ( z , r ) -coordinates

ρ s t 2 u z ( s ) = μ s ( r 1 r ( r r u z ( s ) ) + r 2 ϕ 2 u z ( s ) + z 2 u z ( s ) ) (68)

where u z ( s ) = u z ( s ) ( t , r , ϕ , z ) . In view of (67), oblate spheroidal coordinates are useful. For their definition, we refer to Spencer and Moon [22] . We need to solve the equations R 0 2 = a 2 cosh 2 η 0 and R 0 2 / 2 = a 2 sinh 2 η 0 . This results in a 2 = R 0 2 ( 1 2 ) < R 0 2 and 0 < η 0 = arcoth ( ) . We then let

η ( η 0 δ η / 2 , η 0 + δ η / 2 )

model the finite thickness of the membrane. This is needed because the Helmholtz operator cannot be restricted to the η = η 0 iso-surface in oblate spheroidal coordinates so that a solvable eigenvalue problem is maintained. This is in contrast to spherical coordinates where we can restrict them to the r = R iso-surface.

The coordinate transform we aim at is established as r = a cosh η sin θ and z = a sinh η cos θ . The result is a Laplacian in oblate spheroidal coordinates [22] . Let Φ be a suitably function depending on the oblate spheroidal variables ( η , θ , ϕ ) . The Laplacian Δ Φ is then given by the expression

Δ Φ = 1 a 2 ( cosh 2 η sin 2 θ ) ( 1 cosh η η ( cosh η Φ η ) + 1 sin θ θ ( sin θ Φ θ ) ) + 1 a 2 ( cosh 2 η sin 2 θ ) 2 Φ ϕ 2 . (69)

It is convenient to define k s 2 ρ s ω 2 / μ s for a stimulus of frequency ω . Quasi-stationarity permits the assumption that the shell oscillates with the same frequency. The resulting Helmholtz equation Δ u s ( z ) + k s 2 u s ( z ) = 0 will be separated in oblate spheroidal coordinates in the next paragraph.

Separation of the Helmholtz equation in oblate spheroidal coordinates

As a separation Ansatz, we take u z ( s ) = O ( η ) P ( ϕ ) A ( θ ) with as yet to-be-determined functions O, P & A. Along the lines indicated elsewhere ( [22] , p. 36), the following equations result,

1 cosh η d d η ( cosh η d O d η ) + ( k s 2 a 2 cosh 2 η p ( p + 1 ) + q 2 cosh 2 η ) O ( η ) = 0, (70)

1 sin θ d d θ ( sin θ d A d θ ) + ( k s 2 a 2 sin 2 θ + p ( p + 1 ) q 2 sin 2 θ ) A ( θ ) = 0 , (71)

d 2 P ( ϕ ) d ϕ 2 + q 2 P ( ϕ ) = 0 (72)

where p > 0 and q > 0 are separation parameters. They will be fixed during the discussion to follow. The thickness ( d s ) of the elastic shell in η-direction enters by noting that η ( η 0 δ η / 2 , η 0 + δ η / 2 ) has a small interval length. We force O ( η ) = O ( η 0 ) 0 to be a non-zero constant in the interval. This corresponds to requiring that the whole shell moves along its thickness in η-direction. At θ = π / 2 , the η-average is consistent with the radial average performed to obtain Equation (64) and Equation (63). Physically, the average under consideration formalizes the requirement that only constant elastic deformations with respect to η occur. The equation (70) then gives rise to an equation for the separation parameter p,

k s 2 a 2 cosh 2 η p ( p + 1 ) + q 2 cosh 2 η = 0. (73)

The above equation can be solved for p and the correct sign in front of the square root of the discriminant is fixed by the requirement p > 0 . The result for p is then

p = p ( q , η ) = 1 + 1 + 4 ( k s 2 a 2 cosh 4 η + q 2 cosh 2 η ) 2 . (74)

The function p depends on η and q in a well-behaved way. The dependence upon η can actually be eliminated. The interval for η has a length δ η suitably small so that η η 0 introduces a negligible error. The dependence upon q is fixed by investigating the polar part P = P ( ϕ ) . The clamping of extracolumella needs to be accounted for as the piston mode of the shell plus extracolumella system has already been massaged into (64). The vibrations on the shell’s surface follow from the compressional coupling conditions defined before in (64) & (63), where the shell had been excluded as immobile element. The clamping at ϕ { β ,2 π β } requires Dirichlet boundary conditions for (72):

P ( ϕ = β ) = 0 = P ( ϕ = 2 π β ) ,

which settles as a pair of L 2 ( [ β ,2 π β ] ) -normalized solutions

Φ q ( ϕ ) = sin ( q ( m ) ( ϕ β ) ) π β , (75)

q ( m ) = m π 2 π 2 β ( m ) . (76)

The q-dependence of the parameter p in (74) can be specified explicitly,

p ( m ) = 1 + 1 + 4 ( ( 2 π 2 β ) 2 k s 2 a 2 cosh 4 η 0 + m 2 π 2 ( 2 π 2 β ) 2 cosh 2 η 0 ) 2 . (77)

The obstruction for further analytic investigation is rooted in the solution of (71). Since the separation parameter p is constant with respect to θ , the azimuthal dependence of u z ( s ) is contained in functions of the form

A m ( θ ) = A m ps p ( m ) q ( m ) ( i k s a , cos θ ) + B m qs p ( m ) q ( m ) ( i k s a , cos θ ) . (78)

The objects denoted by ps p ( m ) q ( m ) ( i λ , x ) and qs p ( m ) q ( m ) ( i λ , x ) are known as (azimuthal-) angular part of oblate spheroidal functions. No closed-form analytic representation is known [24] . Following the notation of Falloon, Abbott, and Wang [25] and taking z , μ , γ , they are solutions of the following ordinary differential equation of Sturm-Liouville type,

d d z ( ( 1 z 2 ) d f d z ) + ( λ ν μ + γ 2 ( 1 z 2 ) μ 2 1 z 2 ) = 0 (79)

We refer to the literature [25] [26] [27] for a further treatment of spheroidal wave functions for non-integer parameters p and q. Here we restrict ourselves to noting that the general solution u z ( s ) for the shell equation can be expressed as

u z ( s ) ( t , z , ϕ ) = m = 1 A m ps p ( m ) q ( m ) ( i ω ρ s μ s R 0 1 2 , z R 0 ) sin ( q ( m ) ( ϕ β ) ) 2 π 2 β exp ( i ω t ) + m = 1 B m qs p ( m ) q ( m ) ( i ω ρ s μ s R 0 1 2 , z R 0 ) sin ( q ( m ) ( ϕ β ) ) 2 π 2 β exp ( i ω t ) (80)

The quantity is defined in (66), q ( m ) is specified in (76) and p ( m ) in (77) where η 0 = arcoth ( 2 + λ s / μ s ) is to be used in these equations. The two constants A m and B m for a polar mode of index m are theoretically fixed by the boundary conditions (64) and (63) as well as by the requirement that the top of the shell does not vibrate after exclusion of the piston mode: u ( t , z = R 0 , ϕ ) = 0 . As of 2017 [28] , no stable computer algebra system seems to be available. As such, more than a specification of the general expression (80) is out of reach.

Elastic ring

Expressing the Laplacian Δ in (56) in cylindrical coordinates, a restriction to r = R 0 , u z ( r ) is the solution of the differential equation

ρ r t 2 u z ( r ) = ( λ r + 2 μ r ) z 2 u z ( r ) + μ r R 0 2 ϕ 2 u z ( r ) . (81)

Quasi-stationarity is employed to factor out an exp ( i ω t ) and thus convert the above wave equation to a Helmholtz equation. The method to re-scale the variable z that has been employed to convert (55) to (68) applies to (81) as well, with the obvious substitutions of Lamé constants. Let us define k r 2 = ω 2 ρ r / μ r and redefine z = z where

2 + λ s μ s > 1. (82)

At z = H , we demand that the ring is clamped. The elastic ring is supposed to terminate by being “glued” into the bony structure of the animal’s skull. In the variable z , the clamping condition is imposed at z = H 0 ( ) 1 H < H . Because of the ring’s axial symmetry, periodic boundary conditions are appropriate for the cylindrical angular variable. Then the u z ( r ) solving (81) is of the form

u z ( r ) ( t , ϕ , z ) = m = C m sin ( k s 2 m 2 R 0 2 ( z + H ) ) e i m ϕ e i ω t . (83)

For m , C m is a constant that can be found by insertion of (80) and (83) into the compressional boundary conditions (63) and (64). Together with the condition u z ( s ) ( t , z = R 0 , ϕ ) = 0 , the constants in { ( A m , B m ) | m } for the shell’s displacement as given in (80) and { C m | m } for the ring’s displacement as given in (83) are the solutions to the inhomogeneous system of linear equations that would be defined by the boundary conditions―but:

Analytic investigations of the resulting equations are out of reach because of the lack of closed form expressions for the spheroidal wave functions ps and qs. Because of serious concerns over the reliability of numerical predictions even for integer-indexed spheroidal wave functions, a numerical treatment is ruled out. So we face the question: How can we still make a bit of progress?

Weak coupling assumption

We recall that amplitudes of the piston vibration of the shell have already been shown to be small compared to the vibration amplitude of the ring; cf. Equation (65). Furthermore, biological tissue consists, to a large fraction, of water. Under small external pressures, water is practically incompressible. Also the generation of excess-pressure contributions to the compressional boundary conditions (63) and (64) should be small. In view of the amplification factor a , we therefore assume that the compressional moduli K r , K s of ring and shell are sufficiently small so as to permit iteration of the boundary conditions (63) and (64), a weak coupling assumption.

That is, we assume | K d s z 2 u z ( r ) | , | K d s z 2 u z ( s ) | a p ( t ) . The assumption corresponds to requiring that the elastic coupling in the two-spring model impacts the overall coupled vibration of shell and ring sub-dominantly as compared to the stimulus’ impact. From a biological viewpoint, this behavior is reasonable because the tympanic structure shall respond predominantly only to the external stimulus p that an animal needs to detect. Furthermore, the similarity in material composition between water and biological tissue hints at the smallness of compressional effects compared in driving the membrane system as compared to the influence of the external pressure p. The boundary conditions (64) and (63) are the understood as differential equations defined on the joint boundary at z = 0 and finally subjected to a Picard-iteration. The iterated analogues of (64) and (63) read

σ s t 2 [ u z ( s ) ] ( k + 1 ) ( t , 0 , ϕ ) = ( K s d s ) z 2 [ u z ( s ) ] ( k ) ( t , 0 , ϕ ) + ( K r d s ) z 2 [ u z ( r ) ] ( k ) ( t , 0 , ϕ ) , (84)

σ r t 2 [ u z ( r ) ] ( k + 1 ) ( t , 0 , ϕ ) = ( K r d s ) z 2 [ u z ( r ) ] ( k ) ( t , 0 , ϕ ) + ( K s d s ) z 2 [ u z ( s ) ] ( k ) ( t , 0 , ϕ ) + a p ( t ) (85)

at z = 0 .

The iteration index k is a non-negative integer. As a starting value for the boundary conditions, we specialize to [ u z ( s ) ] ( 1 ) = 0 = [ u z ( s ) ] ( 1 ) . The fact that the external stimulus is the physical cause for the displacement of the tympanic structure underlies this choice of starting values.

For the 0-th iterate in (84), we find σ s t 2 [ u z ( s ) ] ( 0 ) ( t , ϕ , z = 0 ) = 0 for the shell, which is fulfilled, if A m = B m = 0 for all m so that u s ( z ) ( t , ϕ , z ) = 0 . The result means that apart from the piston mode, which has already been translated into a contribution to the ring’s displacement, the shell exhibits no vibrations. For the ring, (85) produces

σ s t 2 [ u z ( r ) ] ( 0 ) ( t , ϕ , z = 0 ) = a p ( t ) (86)

as 0-th iterate boundary condition. In the quasi-stationary state

[ u z ( r ) ] ( 0 ) exp ( i ω t ) ,

the above equation reduces to inhomogeneous Dirichlet boundary conditions, viz.,

[ u z ( r ) ] ( 0 ) ( t , ϕ , z = 0 ) = a p 0 ρ r d s ω 2 exp ( i ω t ) . (87)

Since the homogeneous part of the boundary conditions in 0-th iteration is Dirichlet, the eigenfrequencies of the ring’s vibration pattern in 0-th order can be obtained, following the standard procedure to solve Dirichlet eigenvalue problems, from the requirement k r 2 m 2 / R 0 2 = π n / H in (83). This is equivalent to requiring that the z-dependent contribution to u z ( r ) in (83) satisfies homogeneous Dirichlet boundary conditions, so that

ω ω n m = μ r ρ r m 2 R 0 2 + λ r + 2 μ r ρ r n 2 π 2 H 2 (88)

where n and m . It remains to solve (56) with 0-th order boundary condition (87) imposed instead of the full (64). The approximate boundary value problem gives rise to approximate solutions to (55) and (56) where, instead of the boundary conditions (64) and (63), the iterated boundary conditions (85) and (84) are used in 0-th order. In doing so, we assume that compressional effects of the membrane tissue is small compared to the external pressure drive.

Due to axial symmetry of the stimulus a p = a p ( t ) , only the m = 0 -mode in (83) contributes. Consequently, the 0-th iterate solution reads

u z ( r ) ( t , z , ϕ ) = R 0 p 0 e i ω t 2 ρ r d s 2 ω 2 sin ( μ r λ r + 2 μ r ρ r ω 2 R 0 2 μ r m 2 μ r R 0 2 ( z + H ) ) sin ( μ r λ r + 2 μ r ρ r ω 2 R 0 2 μ r m 2 μ r R 0 2 H ) , (89)

where a = R 0 / 2 d s has been reinstalled.

Implications for the iteration parameter of the two-spring model

The absence of tools to generate reliable numerical predictions for the coupled motion of shell and ring led to the creation of the two-spring model as an effective theory in which the compressional coupling has been stored in the fit constant k 0 . One possible approach to obtain insights in the strength is by noting k 0 / k r K r / ( λ r + 2 μ r ) in view of the solution (89). Under the additional assumption that the shell vibration supports an analogous sinusoidal behavior as in (89), an analogous estimate results: k 0 / k s K s / ( λ s + 2 μ s ) . In terms of elastic parameters

k 0 k r E r ( 1 + ν r ) 3 E r ( 1 ν r ) = 1 3 1 + ν r 1 ν r , (90)

k 0 k s E s ( 1 + ν s ) 3 E s ( 1 ν s ) = 1 3 1 + ν s 1 ν s . (91)

Since for most materials 0 ν r , ν s 0.5 , the assumption k 0 / k r , k 0 / k s 1 seems quite reasonable.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] van Hemmen, J.L., Christensen-Dalsgaard, J., Carr, C.E. and Narins, P. (2016) Animals and ICE: Meaning, Origin and Diversity. Biological Cybernetics, 110, 237-246.
https://doi.org/10.1007/s00422-016-0702-x
[2] Schnupp, J. and Carr, C.E. (2009) One Hearing with More than One Ear: Lessons from Evolution. Nature Neuroscience, 12, 692-697.
https://doi.org/10.1038/nn.2325
[3] Vedumurdi, A.P. (2018) General Aspects of Sound Localization through Intenally Coupled Ears. PhD Thesis, Technische Universität München, München.
[4] Vossen, C. (2010) Auditory Information Processing in Systems with Internally Coupled Ears. PhD thesis, Technische Universität München, München.
[5] Vedurmudi, A.P., Goulet, J., Christensen-Dalsgaard, J., Young, B.A., Williams, R. and van Hemmen, J.L. (2016) How Internally Coupled Ears Generate Temporal and Amplitude Cues for Sound Localization. Physical Review Letters, 116, 028101.
https://doi.org/10.1103/PhysRevLett.116.028101
[6] C. Vossen, J. Christensen-Dalsgaard, and J. L van Hemmen. Analytical model of internally coupled ears. Journal of the Acoustical Society of America, 128(2):909-918, 2010.
https://doi.org/10.1121/1.3455853
[7] Howe, M.S. (2007) Hydrodynamics and Sound. Cambridge University Press, Cambridge, UK.
[8] Howe, M.S. (2014) Acoustics and Aerodynamic Sound. Cambridge University Press, Cambridge, UK.
[9] Autrum, H. (1940) über Lautäusserungen und Schallwahrnehmungen bei Anthropoden II. Das Richtungshören von Locusta und Versuch einer Hörtheorie für Tympanalorgane vom Locustidentyp. Zeitschrift für Vergleichende Physiologie, 28, 326-352.
[10] Autrum, H. (1942) Schallempfang bei Tier und Mensch. Naturwissenschaften, 30, 69-85.
https://doi.org/10.1007/BF01475622
[11] Vedurmudi, A.P., Christensen-Dalsgaard, J. and van Hemmen, J.L. (2018) Modeling Underwater Hearing and Sound Localization in the Aquatic Frog Xenopus laevis. Journal of the Acoustical Society of America, Submitted.
[12] Vedurmudi, A.P., Young, B.A. and van Hemmen, J.L. (2016) Internally Coupled Ears: Mathematical Structures and Mechanisms Underlying ICE. Biological Cybernetics, 110, 237-246.
https://doi.org/10.1007/s00422-016-0696-4
[13] Fletcher, N.H. (1992) Acoustic Systems in Biology. Oxford University Press, Oxford, UK.
[14] Saunders, J.C., Duncan, R.K., Doan, D.E. and Werner, Y.L. (2000) The Middle Ear of Reptiles and Birds. In: Dooling, R.J., Fay, R.R. and Popper, A.N., Eds., Comparative Hearing: Birds and Reptiles, Springer, New York, Chapter 2.
https://doi.org/10.1007/978-1-4612-1182-2_2
[15] Timoshenko, S.P. and Goodier, J. (1969) Theory of Elasticity. 3rd Edition, McGraw-Hill, New York.
[16] Timoshenko, S.P. and Woinowsky-Krieger, S. (1959) Theory of Plates and Shells. 2nd Edition, McGraw-Hill, New York.
[17] Van Hemmen, J.L. and Leibold, C. (2007) Elementary Excitations of Biomembranes: Differential Geometry of Undulations in Elastic Surfaces. Physics Reports, 444, 51-99.
https://doi.org/10.1016/j.physrep.2006.12.007
[18] Hassani, S. (2013) Mathematical Physics: An Introduction to Its Foundations. 2nd Edition, Springer, Berlin.
https://doi.org/10.1007/978-3-319-01195-0
[19] Van der Waerden, B.L. (1993) Algebra I. 9th Edition, Springer, Berlin.
https://doi.org/10.1007/978-3-642-85527-6
[20] Van der Waerden, B.L. (1993) Algebra II. 6th Edition, Springer, Berlin.
https://doi.org/10.1007/978-3-642-58038-3
[21] Colbert, E. (1946) The Eustachian Tubes in the Crocodilia. Copeia, 12, 12-14.
https://doi.org/10.2307/1438813
[22] Spencer, P. and Moon, D. (1971) Field Theory Handbook: Including Coordinate Systems, Differential Equations and their Solutions. Springer, Berlin.
[23] Lautrup, B. (2004) Physics of Continuous Matter: Exotic and Everyday Phenomena in the Macroscopic World. IoP Publishing, Bristol.
[24] Skudrzyk, E. (1971) Foundations of Acoustics: Basic Mathematics and Basic Acoustics. Springer, Berlin.
https://doi.org/10.1007/978-3-7091-8255-0
[25] Falloon, P., Abbott, P. and Wang, J. (2002) Theory and Computation of the Spheroidal Wave Functions.
https://arxiv.org/ftp/math-ph/papers/0212/0212051.pdf
[26] Meixner, J. and Schäfke, R. (1954) Mathieusche Funktionen und Sphäroidfunktionen. Springer, Berlin.
https://doi.org/10.1007/978-3-662-00941-3
[27] Meixner, J., Schäfke, R. and Wolf, G. (1980) Mathieu Functions and Spheroidal Functions and Their Mathematical Foundations (Further Studies). Lecture Notes in Mathematics 837. Springer, Berlin.
https://doi.org/10.1007/BFb0096194
[28] Zhao, L. (2017) Spherical and Spheroidal Harmonics: Examples and Computations. Master’s Thesis, Ohio State University, Columbus.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.