Mimicking the active cochlea with a fluid-coupled array of subwavelength Hopf resonators

TL;DR

This paper introduces a fluid-coupled array of subwavelength Hopf resonators that mimics the active cochlea's frequency separation and amplification, advancing acoustic metamaterials for auditory applications.

Contribution

It presents a novel design of an acoustic metamaterial with non-linear amplification, modeling cochlear behavior using a size-graded array of resonators with modal decomposition.

Findings

Successfully reproduces cochlear frequency separation

Implements non-linear amplification similar to cochlear amplifier

Provides a physically-based modal analysis of the resonator array

Abstract

We present a design for an acoustic metamaterial that mimics the behaviour of the active cochlea. This material is composed of a size-graded array of cylindrical subwavelength resonators, has similar dimensions to the cochlea and is able to reproduce the frequency separation of audible frequencies. Non-linear amplification is introduced to the model in order to replicate the behaviour of the cochlear amplifier. This formulation takes the form of a fluid-coupled array of Hopf resonators. We seek solutions based on a modal decomposition, so as to retain the physically-derived coupling between resonators.

Equations

$$\frac{\mathrm{d}z}{\mathrm{d}t}=(\mu +i{\omega }_0)z-|z{|}^{2}z+F,$$

$$v_0=\sqrt{\frac{{\kappa }_0}{{\rho }_0}},v=\sqrt{\frac{\kappa }{\rho }}.$$

$$\delta :=\frac{\rho }{{\rho }_0}.$$

$$p(x,t)=\mathrm{Re}\left(\sum _n{\alpha }_n(t)u_n(x)\right).$$

$$\begin{cases}\left(\Delta+\tfrac{\omega^{2}}{v_{0}^{2}}\right)u(x,\omega)=0,&\text{for }(x,\omega)\in\mathbb{R}^{2}\setminus\overline{D}\times\mathbb{C},\\ \left(\Delta+\tfrac{\omega^{2}}{v^{2}}\right)u(x,\omega)=0,&\text{for }(x,\omega)\in D\times\mathbb{C},\\ u_{+}-u_{-}=0,&\text{for }(x,\omega)\in\partial D\times\mathbb{C},\\ \delta\frac{\partial u}{\partial\nu}\big{|}_{+}-\frac{\partial u}{\partial\nu}\big{|}_{-}=0,&\text{for }(x,\omega)\in\partial D\times\mathbb{C},\end{cases}$$

$$\underset{|x|\to \infty }{\lim }|x{|}^{1/2}\left(\frac{\partial }{\partial |x|}-i\frac{\omega }{v_0}\right)u(x,\omega )=0.$$

$$u(x,\omega )=\{\begin{array}{ll}{\mathcal{S}}_D^{\omega /v_0}[\psi ](x),& (x,\omega )\in {\mathbb{R}}^2\setminus \overline{D}\times \mathbb{C},\\ {\mathcal{S}}_D^{\omega /v}[\varphi ](x),& (x,\omega )\in D\times \mathbb{C},\end{array}$$

$${\mathcal{S}}_D^{\omega }[\phi ](x):={\int }_{\partial D}{\Gamma }^{\omega }(x-y)\phi (y)d\sigma (y),x\in \partial D,\phi \in L^2(\partial D),\omega \in \mathbb{C},$$

$$\mathcal{N}[\phi ]:=\mu \phi -\beta |\phi {|}^2\phi .$$

$$p(x,t)=\mathrm{Re}\left(\sum _{n=1}^N{\alpha }_n(t)u_n(x)\right),$$

$$f(t){\mathcal{X}}_Q(x)\simeq f(t)\sum _{n=1}^N{\mathcal{F}}_n{u}_n(x),$$

$$\sum_{n=1}^{N}\big{(}\alpha_{n}^{\prime\prime}(t)+\omega_{n}^{2}\alpha_{n}(t)\big{)}u_{n}(x)=f(t)\sum_{n=1}^{N}\mathcal{F}_{n}u_{n}(x)+\mathcal{N}\left[\sum_{n=1}^{N}\alpha_{n}^{\prime}(t)u_{n}(x)\right]\mathcal{X}_{D}(x).$$

$${\mathbf{\alpha }}^{\prime \prime }+\mathbf{\Lambda }\mathbf{\alpha }-\mu {\mathbf{\alpha }}^{\prime }+\beta \mathbf{N}({\mathbf{\alpha }}^{\prime })=f(t)\mathbf{F},$$

$$\mathbf{N}(\mathbf{z}):={\gamma }^{-1}{\left[{\text{bigintsss}}_D{\left|\sum _{n=1}^N{z}_n{u}_n(x)\right|}^2\left(\sum _{n=1}^N{z}_n{u}_n(x)\right)\overline{u_j(x)}\mathrm{d}x\right]}_{j=1,...,N}.$$

$${\alpha }_k^{\prime \prime }+{\omega }_k^2{\alpha }_k-\mu {\alpha }_k^{\prime }+\hat{\beta }|{\alpha }_k^{\prime }{|}^2{\alpha }_k^{\prime }=f(t){\mathcal{F}}_k,$$

$$\begin{array}{c}{\Omega }_k^c:=\sqrt{\mathrm{Re}({\omega }_k{)}^2-\mathrm{Im}({\omega }_k{)}^2},\\ {\mu }_k^c:=\frac{-2\mathrm{Re}({\omega }_k)\mathrm{Im}({\omega }_k)}{{\Omega }_k^c},R_k^c=\sqrt{\frac{\mu -{\mu }_k^c}{\hat{\beta }}}\frac{1}{{\Omega }_k^c}.\end{array}$$

$$J=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -({\Omega }_k^c{)}^2& {\mu }_k^c{\Omega }_k^c& \mu & 0\\ -{\mu }_k^c{\Omega }_k^c& -({\Omega }_k^c{)}^2& 0& \mu \end{array}\right],$$

$$\lambda =\frac{1}{2}\left(\mu \pm \sqrt{-4({\Omega }_k^c{)}^2\pm 4i{\Omega }_k^c{\mu }_k^c+{\mu }^2}\right).$$

$$(-2i{\Omega }_k^c{R}_k^{\prime }-({\Omega }_k^c{)}^2{R}_k+{\omega }_k^2{R}_k)-\mu (R_k^{\prime }-i{\Omega }_k^c{R}_k)+\hat{\beta }(({\Omega }_k^c{)}^2{R}_k^2{R}_k^{\prime }-i({\Omega }_k^c{)}^3{R}_k^3)=0,$$

$${\alpha }_k(t)=R_k(t)e^{-i{\Omega }_k^{c}t+i{\psi }_k(t)},{\alpha }_k^{\prime }(t)=-i{\Omega }_k^c{R}_k(t)e^{-i{\Omega }_k^{c}t+i{\psi }_k(t)}.$$

$$\left(R_k^{\prime }+i{\psi }_k^{\prime }{R}_k\right)e^{-i\Omega t+i{\psi }_k}=0.$$

$$\left(-i{\Omega }_k^c{R}_k^{\prime }+{\Omega }_k^c{\psi }_k^{\prime }{R}_k\right)e^{-i\Omega t+i{\psi }_k}=ϵi\left(-{\Omega }_k^c{R}_k+({\Omega }_k^c{)}^3{R}_k^3\right)e^{-i\Omega t+i{\psi }_k}.$$

$$R_k^{\prime }$$

$${\psi }_k^{\prime }$$

$$R_k$$

$${\psi }_k$$

$${\tilde{R}}_k^{\prime }$$

$${\tilde{\psi }}_k^{\prime }$$

$$\begin{array}{rl}{h}_1({\tilde{R}}_k,{\tilde{\psi }}_k,t)& ={\int }_0^t\left({\tilde{R}}_k-({\Omega }_k^c{)}^2{\tilde{R}}_k^3\right){\sin }^2({\Omega }_k^{c}t-{\tilde{\psi }}_k)\mathrm{d}t\\ & -\left[\frac{{\Omega }_k^c}{2\pi }{\int }_0^{2\pi /{\Omega }_k^c}\left({\tilde{R}}_k-({\Omega }_k^c{)}^2{\tilde{R}}_k^3\right){\sin }^2({\Omega }_k^{c}t-{\tilde{\psi }}_k)\mathrm{d}t\right]t.\end{array}$$

$${\tilde{R}}_k^{\prime }=\frac{1}{2}ϵ\left({\tilde{R}}_k-({\Omega }_k^c{)}^2{\tilde{R}}_k^3\right),{\tilde{\psi }}_k^{\prime }=0.$$

$$R_k(t)=\frac{1}{\sqrt{R_k(0{)}^{-2}{e}^{-ϵt}+({\Omega }_k^c{)}^2\left(1-e^{-ϵt}\right)}}+O(ϵ).$$

$${\left(({\Omega }_k^c{)}^2-{\Omega }^2\right)}^2{R}_k^2+{\left(-{\mu }_k^c{\Omega }_k^c{R}_k+\mu \Omega R_k-\hat{\beta }{\Omega }^3{R}_k^3\right)}^2=F^2|{\mathcal{F}}_k{|}^2.$$

Full text

Mathematics subject classification: 35Q92, 37N25, 35J05

Keywords: subwavelength resonance, acoustic metamaterials, cochlear mimicry, coupled Hopf resonators, active cochlear mechanics, non-linear amplification

1 Introduction

1.1 The cochlea

The cochlea is the part of the ear responsible for transforming physical vibrations into neural signals. It does so with remarkable efficacy, capturing a large range of volumes with impressive frequency resolution. At its simplest, the cochlea is a long tube that is filled with fluid and divided in two by the cochlear partition [1]. This partition is set in motion by acoustic waves in the cochlear fluid and contains receptor cells which respond to this motion [2, 3, 4]. These cells, known as hair cells, can be split into two categories, referred to as inner and outer due to their positioning. The former are responsible for converting motion of the partition into electrical signals, while the latter contain a motor protein and are involved in the cochlea’s non-linear amplification mechanism.

The cochlea’s mechanism for resolving frequency is based on the fact that several of its physical parameters (such as the stiffness and width of the basilar membrane, upon which the cochlear partition is mounted) are graded along its length. This property means that different tones will lead to maximal excitation at different points [5]. High frequency sounds are detected at the base of the cochlea while sounds of lower frequencies are detected closer to the apex. This frequency separation mechanism is remarkably effective, meaning that even non-musicians can detect differences of less than 0.5% [6, 7].

Humans are also able to comfortably hear sounds with amplitudes that range over an impressive six orders of magnitude. This is thanks to the cochlea’s non-linear amplification response, known as the cochlear amplifier. This is such that a compressive non-linearity occurs, whereby quieter sounds are amplified much more greatly than louder ones. Understanding the details of the mechanism responsible for the cochlear amplifier represents the most significant open problem in the study of hearing. Recently, the most prevalent models have been based on a critically-poised Hopf resonator, the details of which are discussed in Section 1.3.

1.2 An acoustic metamaterial

Metamaterials are broadly defined as synthetic structures which have remarkable properties that do not generally exist in conventional materials. Typically composed of a repeating small structure, acoustic metamaterials have been successfully used for a variety of wave-control applications [8, 9, 10]. In this work we will study a gradient acoustic metamaterial, where the material parameters (in particular, the size of the small, repeating unit) vary with position. It is known that metamaterials of this kind can perform spatial frequency separation. This phenomenon is often known as rainbow trapping, particularly in the context of optics where it was first proposed by Tsakmakidis et al. [11]. It has since been observed a in variety of different settings including acoustic metamaterials [12, 13, 14, 15], amongst others [16, 17, 18, 19].

In this work, we study an acoustic metamaterial composed of an array of high-contrast subwavelength resonators, graded in size. These are compressible objects with material parameters (density and bulk modulus) that differ greatly from the background medium and, as a result, have a resonant frequency that corresponds to a wavelength much greater than the object’s physical size. The classical example is an air bubble in water (known as a Minnaert bubble [20]) but the theory holds for any pair of materials with a sufficiently great parameter contrast. An asymptotic analysis of a graded array of such resonators was recently performed by Ammari and Davies [13], where it was shown that with a suitably-chosen gradient the structure reproduced the spatial frequency separation exhibited by the cochlea [5].

Rupin et al. [21] have recently produced promising experimental results, demonstrating that a graded array of cylindrical resonators can mimic both the cochlea’s frequency separation and non-linear amplification properties. Their set-up was scaled up for practicality and employed an array of cylindrical quarter-wavelength resonators with a microphone-and-speaker amplification system. In this work we study a structure with similar dimensions to the cochlea, consisting of gas-filled cylinders surrounded by fluid. We take a cross section of this structure so that we study the problem of an acoustic wave scattered in a plane by a size-graded array of circles.

As an aside, a parallel source of inspiration for the gradient metameterial considered in this work is the assortment of cochlear models based on an arrays of harmonic oscillators [22, 23, 24, 25, 26, 27]. The compliance of the cochlear partition is largely accounted for by the basilar membrane which is composed of soft, elastic tissue reinforced by strong collagen fibres across the cochlea’s width (Figure 1) [28, 29, 3]. As a result, if one considers the membrane in sections then the mechanical coupling between them is relatively small and each can be modelled as an individual harmonic oscillator. This motivates discrete cochlear models based on graded arrays of, for example, masses on springs [22] or elastic beams [23]. In fact, many of the earliest theories of hearing were based on a version of this principle [30], such as Helmholtz’ famous “piano strings” analogy [31]. Since Minnaert-type resonators behave as harmonic oscillators [32, 20], the structure studied here can be viewed as a manifestation of these models.

The choice of resonators made here is significant since it allows for the design of a structure that has similar dimensions to a cochlea and responds to audible frequencies. Throughout this paper we will run simulations on an array of 22 resonators, arranged in a linear configuration that measures $32$\mathrm{mm} in length and $0.28$\mathrm{mm} at its widest. We use the material parameters of air and water for the resonators and surrounding fluid, respectively, in order to demonstrate that the material properties (in particular, the contrast) required to give the desired subwavelength resonant response are exhibited by everyday substances. Bubbly structures of this kind have previously been constructed e.g. by injecting air bubbles into silicone-based polymers [33, 34].

Modelling the interactions between the array of resonators, which occur on subwavelength scales, is a challenging problem. We will use a modal-decomposition approach to study the wave-propagation problem [35, 36]. The structure’s spatial eigenmodes will be found using layer-potential formulations [37] before these profiles are used as a basis to study the behaviour of the system when appropriate non-linear amplification is introduced to the resonators.

1.3 Hopf resonators in cochlear mechanics

While there is still much to be learnt about the cochlear amplifier, researchers are guided by a handful of valuable observations. It is understood, firstly, that this mechanism acts as a negative-damping mechanism [38, 39, 40]. Secondly, it has been observed that non-linear amplification occurs with a one-third power law, suggesting the presence of a cubic non-linearity [6].

Hopf resonators have become popular objects to study in the field of non-linear cochlear mechanics thanks to their remarkable ability to account for the key properties that typify cochlear behaviour [41, 6, 42, 43, 40, 44, 25, 24, 45]. The normal form of a single Hopf resonator $z=z(t):\mathbb{R}\to\mathbb{C}$ in the complex plane is given by the forced differential equation

[TABLE]

where $F=F(t)$ is the forcing term and $\omega_{0}$ and $\mu$ are real parameters. This system is a resonator in the sense that the absolute value of the response $z$ is greatest when the forcing $F$ occurs with frequency $\omega_{0}$. In cochlear models, $z$ is some variable which characterises the system’s state.

The parameter $\mu$ is the bifurcation parameter. For $\mu<0$ the unforced system ($F=0$) has a stable equilibrium at $z=0$ whereas when $\mu>0$ this equilibrium is unstable and there exists a stable limit cycle given by $z(t)=\sqrt{\mu}e^{i\omega_{0}t}$. This birth of a limit cycle is typical of a (supercritical) Hopf bifurcation, which is formally characterised by a conjugate-pair of linearised eigenvalues crossing the imaginary axis. Writing the unforced system (1.1) in terms of its real and imaginary parts and linearising about the fixed point at zero gives a system whose Jacobian matrix has eigenvalues $\lambda=\mu\pm i\omega_{0}$. These eigenvalues clearly cross the imaginary axis when $\mu$ passes zero. For further details see e.g. [46, 47, 35].

The greatly enhanced response for frequencies close to $\omega_{0}$ is able to account for the cochlea’s frequency selectivity. The cubic non-linearity in (1.1) is able to reproduce the one-third power law of the cochlea: when $\mu$ is small, and the system is close to bifurcation, we have that $|z|\approx|F|^{1/3}$ for frequencies close to resonance.

One of the earliest pieces of evidence supporting the active nature of the cochlea was the observation that the ear emits sounds (known at otoacoustic emissions) as part of its response [48, 49]. The existence of stable limit cycles, for certain parameter values, predicts this behaviour [43, 44, 50].

A further symptom of the non-linearity that exists in the cochlea is the behaviour that is observed under the influence of a signal composed of two distinct tones. Firstly, when the ear is excited by such a stimulus two-tone suppression occurs. That is, the frequency spectrum of the response contains the expected two amplitude peaks, however, these are smaller than each would be in the absence of the other tone [51]. Further, in this situation the ear also detects additional tones, variously known as combination tones, distortion products or Tartini’s tones [31, 52, 53]. Close to bifurcation, the non-linearity in (1.1) gives products that can account for these phenomena [53, 44].

In this work, we will introduce a Hopf-type non-linearity directly to the wave-propagation problem by supposing that the resonators are equipped with an appropriate forcing mechanism (as was realised by [21]). Based on an eigenmode decomposition, we will explore the Hopf-type behaviour of this system and show that the crucial cochlea-like properties of (1.1) are retained by the coupled subwavelength structure.

2 Scattering by coupled subwavelength resonators

It is known that the frequency-location (tonotopic) map in the cochlea is exponential [5] and it was shown in [13] that if an array of subwavelength resonators has a similar size-grading, we reproduce this phenomenon. Thus, we will consider a domain $D$ in $\mathbb{R}^{2}$ which is the disjoint union of $N\in\mathbb{N}$ circular subdomains $\{D_{1},\ldots,D_{N}\}$ with each successive radius defined as $R_{n+1}=sR_{n}$, for some $s>1$.

We denote by $\rho$ and $\kappa$ the density and bulk modulus of the interior of the resonators, respectively, and use $\rho_{0}$ and $\kappa_{0}$ for the corresponding parameters for the background fluid (which occupies $\mathbb{R}^{2}\setminus\overline{D}$). We may then denote the acoustic wave speeds in $\mathbb{R}^{2}\setminus\overline{D}$ and $D$, respectively, by

[TABLE]

We also introduce the dimensionless contrast parameter

[TABLE]

Since we want our structure to be of similar dimensions to the cochlea and exhibit a resonant response to audible frequencies, we need subwavelength resonant modes to exist. This is known to occur in the case that $\delta\ll 1$ [20, 54, 13]. For the simulations performed here we use material parameters corresponding to air-filled resonators surrounded by water, giving that $\delta\approx 10^{-3}$.

We will use a modal decomposition to analyse the wave propagation within the structure [35, 36]. That is, we wish to express the acoustic pressure $p=p(x,t)$ at position $x$ and time $t$ in the form

[TABLE]

Separating variables in the unforced, linear wave equation yields a spatial Helmholtz problem given by

[TABLE]

where $\omega\in\mathbb{C}$ is the separation constant, $\frac{\partial}{\partial\nu}$ denotes the outward normal derivative and the subscripts + and - are used to denote evaluation from outside and inside $\partial D$, respectively. We must also insist that $u(\cdot,\omega)$ satisfies the Sommerfeld radiation condition

[TABLE]

This condition is required to ensure that the solution represents outgoing waves (rather than incoming from infinity) and gives the well-posedness of (2.4). Since we study the problem in unbounded space, energy is lost to the far-field.

In light of the fact that (2.4) contains no forcing or amplification, we define a resonant frequency and associated eigenmode (or resonant mode) to be solutions $(\omega,u(\cdot,\omega))\in\mathbb{C}\times H_{loc}^{1}(\mathbb{R}^{2})$ of (2.4). Here, $H_{loc}^{1}(\mathbb{R}^{2})$ is the space of functions that, on every compact subset of $\mathbb{R}^{2}$, are square integrable and have a weak first derivative that is also square integrable. We are interested in solutions where $\omega$ is small and the resonators are much smaller than the wavelength of the associated radiation, such solutions are referred to as subwavelength modes.

2.1 Layer-potential approach

We solve (2.4) using a layer-potential approach [37]. This is based on representing the solution, in terms of some surface potentials $\phi,\psi\in L^{2}(\partial D)$, as

[TABLE]

where $\mathcal{S}_{D}^{\omega}$ is the Helmholtz single layer potential associated with the domain $D$. This integral operator is defined as

[TABLE]

where $\Gamma^{\omega}$ is the outgoing (i.e. satisfying the Sommerfeld radiation condition) fundamental solution to the Helmholtz operator $\Delta+\omega^{2}$ in $\mathbb{R}^{2}$ [55]. The value of this approach is that solving (2.4) is reduced to finding $\phi,\psi$ such that the two transmission conditions on $\partial D$ hold [37, 54].

An asymptotic analysis of the resonant frequencies and eigenmodes, based on their layer-potential representations (2.6), was performed in [13]. It was shown that a system of $N$ coupled resonators has $N$ subwavelength resonant modes $u_{1}(x),\dots,u_{N}(x)$ and corresponding resonant frequencies $\omega_{1},\dots,\omega_{N}$ with positive real part.

Using the layer-potential representation (2.6), we are able to find the subwavelength resonant modes numerically by expanding the functions $\phi$ and $\psi$ in terms of Fourier bases on the boundary of each resonator. This is particularly convenient in the case of circular resonators since the evaluation of $\mathcal{S}_{D}^{\omega}$ on Fourier modes has closed-form expressions, the details of which are given in Appendix A.

The $N$ subwavelength eigenmodes take the form of increasingly oscillating coupled patterns (Figure 2). The profiles share several similarities with the response of the basilar membrane. Each mode has a position of maximal amplitude (sometimes known as the resonant place) that depends on its resonant frequency [5, 13]. In the region of the resonant place, the wavelength decreases as the amplitude peaks, before the solution dies away quickly. This was similarly observed by [21] in simulations using point scatterers. The use of resonators with non-zero radii, and the fact that each mode is approximately constant on each resonator [13], mean the profiles look less smooth here.

The design studied here (with $R_{1}=0.1$\mathrm{mm}$$ and $s=1.05$) is chosen so that resonant frequencies (Figure 2) have real parts which fall within the range of audible frequencies (often quoted as 20 Hz - 20 kHz). The negative imaginary parts denote the positive rate of attenuation associated with each mode, due to energy being lost to the far field (note that the assumed time harmonicity in (2.5) is $e^{-i\omega t}$). We will see in Section 3 that the amplification required to illicit a Hopf bifurcation depends on this quantity.

The structure will also have higher order resonant modes at frequencies which correspond to wavelengths similar to the size of the resonators, or bigger. Given the physical dimensions and wavelengths of the problem we are interested in, we focus our attention on the subwavelength modes as these will dominate the behaviour of the system.

3 Non-linear system

We now wish to introduce appropriate non-linear amplification to the model. As discussed in Section 1.3, the canonical form of a Hopf resonator is able to account for the important properties of the cochlear amplifier. This suggests adding amplification based on a non-linearity of the form

[TABLE]

These two terms, respectively, account for the negative damping and cubic non-linearity that we said our amplification should include.

An important consideration, when choosing to introduce amplification, is the stability of the system. For example, Rupin et al. [21] used a formulation whereby amplification closely resembling $\mathcal{N}[p]$ was added. In order for this formulation to be stable, it was necessary to design a set-up that switched off the amplification if the pressure exceeded a threshold value. Conversely, there exist a number of formulations which are stable without this thresholding. Examples include variants of $\mathcal{N}[\partial_{t}p]$ used in the artificial cochlear devices of Joyce and Tarazaga [40, 56, 57] and the $\mu\partial_{t}p-\beta|p|^{2}\partial_{t}p$ term considered by Duke and Jülicher [44].

In this work we will study the system produced by introducing amplification of the form $\mathcal{N}[\partial_{t}p]$ to the resonators. As we shall see in Section 3.2.2, this system is stable without the need to impose a pressure threshold. Furthermore, there is evidence which suggests that hair cell stimulation (by stereocilia displacement) is dependent on membrane velocity [58], suggesting that any amplification should be a function of $\partial_{t}p$.

Since the $N$ subwavelength modes are expected to dominate, we truncate the expansion (2.3), seeking a solution of the form

[TABLE]

for some complex-valued functions of time $\alpha_{1}(t),\dots,\alpha_{N}(t)$. We will also project the forcing term onto the space spanned by the subwavelength eigenmodes. If $Q\subset\mathbb{R}^{2}$ is a compact set on which the forcing is applied (and $\mathcal{X}_{Q}$ is the characteristic function of $Q$) then we decompose the forcing as

[TABLE]

where $\mathcal{F}_{n}:=(\mathcal{X}_{Q},u_{n})_{2,Q}$.

In light of the transmission properties (across $\partial D$) that the eigenmodes inherit from (2.4), we reach the problem

[TABLE]

3.1 Modal system

Our approach to studying (3.4) will be to take the $L^{2}(D)$ product with $u_{m}$ for $m=1,\dots,N$ to reach a coupled system of $N$ ordinary differential equations. Define the matrix $\gamma\in\mathbb{C}^{N\times N}$ as $\gamma_{ij}:=(u_{i},u_{j})_{L^{2}(D)}$. An important property is that, thanks to the linear independence of the eigenmodes, $\gamma$ is invertible [13]. If we write $\alpha_{1},\dots,\alpha_{N}$ in the column vector $\boldsymbol{\alpha}$ then the modal system is described by

[TABLE]

where $\mathbf{\Lambda}\in\mathbb{C}^{N\times N}$ is diagonal with entries $\omega_{n}^{2}$, $\mathbf{F}\in\mathbb{C}^{N}$ is the vector of forcing constants $\mathcal{F}_{n}$ and $\mathbf{N}:\mathbb{C}^{N}\to\mathbb{C}^{N}$ is the non-linear function defined as

[TABLE]

In much of what follows we will be interested in the case when $f(t)$ is harmonic with frequency $\Omega$. In this case, we can approximate the solution to (3.5) using a harmonic balance approach [35, 59, 36, 60]. That is, if $f(t)=Fe^{-i\Omega t}$, for $F,\Omega\in\mathbb{R}$, then we may approximate the steady-state solutions to (3.5) as $\alpha_{k}(t)=X_{k}e^{-i\Omega t+i\psi_{k}}$ for amplitudes $X_{k}\in\mathbb{R}$ and phase delays $\psi_{k}\in\mathbb{R}$. Making this substitution leads to a system of coupled cubic equations that can be solved numerically.

In Figure 3 we see that, as is to be expected, that as the forcing frequency is varied each mode is excited much more greatly in the vicinity of the associated resonant frequency (in spite of the coupling within the non-linearity (3.6)). This motivates an approximate system whereby, if the system is forced at a frequency close to one of the resonant frequencies, we assume that only that mode is excited.

3.2 Single-mode approximation

When the forcing frequency $\Omega$ is close to one of the resonant frequencies $\omega_{k}$ we approximate the solution to (3.4) by assuming that only the corresponding mode $u_{k}$ is excited. In such a regime, we take the $L^{2}(D)$ product of (3.4) with $u_{k}$ to yield the equation

[TABLE]

where $\hat{\beta}:=\beta\|u_{k}\|_{4,D}^{4}/\|u_{k}\|_{2,D}^{2}$.

3.2.1 Hopf bifurcation

At this point, we pause to explore the Hopf-type behaviour that is exhibited by our model. In the case that $f=0$, we see that (3.7) has a periodic solution $\alpha_{k}(t)=R_{k}^{c}e^{-i\Omega_{k}^{c}t}$ provided that $\mu\geq\mu_{k}^{c}$, where

[TABLE]

This birth of a limit cycle is typical of a Hopf bifurcation. A Hopf bifurcation is characterised by a conjugate pair of linearised eigenvalues crossing the imaginary axis [46]. Decomposing $\alpha_{k}$ into its real and imaginary parts, we can write (3.7) as a four-dimensional system of first-order ordinary differential equations. Linearising this system around the fixed point at $\alpha_{k}=0$ gives the Jacobian matrix

[TABLE]

which has eigenvalues given, using the notation of (3.8), by

[TABLE]

When $\mu=\mu_{k}^{c}$, the eigenvalues of $J$ are $\lambda=\pm i\Omega_{k}^{c},\mu_{k}^{c}\pm i\Omega_{k}^{c}$. It can also be shown that $\frac{\>\mathrm{d}}{\>\mathrm{d}\mu}\operatorname{Re}(\lambda)>0$ meaning that the pair of eigenvalues cross the imaginary axis (from left to right) as $\mu$ passes the critical value.

In order to visualise the local stability of this limit cycle and the fixed point at $\alpha_{k}=0$, we allow the radius $R_{k}$ to be a slowly varying function of $t$. That is, we use the ansatz $\alpha_{k}=R_{k}(t)e^{-i\Omega_{k}^{c}t}$ in (3.7) and disregard any terms containing either $R_{k}^{\prime\prime}$ or products of $R_{k}^{\prime}$. This approach leads to the equation

[TABLE]

which is linear in $R_{k}^{\prime}$. The phase planes for $\mu>\mu_{k}^{c}$ and $\mu<\mu_{k}^{c}$ (Figure 4) demonstrate that as $\mu$ passes the critical value a stable limit cycle (described by (3.8)) is born out of the stable equilibrium at the origin, as is typical of a (supercritical) Hopf bifurcation.

3.2.2 Stability of solutions

An important consideration, when choosing an appropriate non-linearity, is the stability of the system. We explore the stability of unforced solutions to (3.7) using a technique known as averaging [35, 61, 60]. The analysis in Section 3.2.1 showed that there is a locally stable limit cycle (when $\mu>\mu_{k}^{c}$) but it is valuable to understand what happens if $R_{k}(0)$ is further away from $R_{k}^{c}$. We begin with the ansatz

[TABLE]

We assume for this analysis that $\hat{\beta}=:\epsilon>0$ and that $\mu=\mu_{k}^{c}+\epsilon$.

Differentiating the first expression of (3.12) and substituting into the second yields

[TABLE]

Further, substituting (3.12) into (3.7) gives

[TABLE]

We may take the real parts of (3.13) and (3.14) and solve for $R_{k}^{\prime}$ and $\psi_{k}^{\prime}$ to give

[TABLE]

We now make a near-identity transformation in order to express $R_{k}$ and $\psi_{k}$ in terms of their average values over the interval $(t-\pi/\Omega_{k}^{c},t+\pi/\Omega_{k}^{c})$, which we denote by $\tilde{R}_{k}$ and $\tilde{\psi}_{k}$. This transformation has the form

[TABLE]

where $h_{1}$ and $h_{2}$ should be chosen in order to simplify the equations for $\tilde{R}_{k}$ and $\tilde{\psi}_{k}$ as much as possible. This substitution leads to the equations

[TABLE]

Ideally, we would like to choose $h_{1}$ so that it cancels with the other $O(\epsilon)$ term in (3.19). However, this antiderivative might grow in time meaning the expansion (3.17) will not be valid for large $t$. Instead, we take $h_{1}$ as the antiderivative minus a linear term that grows with the average value [61], that is

[TABLE]

After substitution of (3.21) into (3.19), we make an approximation in the spirit of the “averaging” methodology [35, 61, 60]. We will assume that the integral in the second term of (3.21) can be well-approximated by taking the value of $\tilde{R}_{k}$ and $\tilde{\psi}_{k}$ as constant over a cycle of oscillation, leaving a simple trigonometric integral.

We choose $h_{2}$ similarly and find that, up to an error of order $O(\epsilon^{2})$,

[TABLE]

Solving by separation of variables gives that

[TABLE]

Crucially, for any $R_{k}(0)>0$ it holds that $R_{k}(t)\to R_{k}^{c}$ as $t\to\infty$, demonstrating that this limit cycle is asymptotically stable.

3.2.3 Pure-tone response

Consider the case of an incoming signal that consists of a single pure tone at frequency $\Omega$, that is, $f(t)=Fe^{-i\Omega t}$ for $F,\Omega\in\mathbb{R}$, where $\Omega$ is close to $\omega_{k}$. Using the harmonic balance ansatz $\alpha_{k}(t)=R_{k}e^{-i\Omega t+i\psi_{k}}$ and finding the complex modulus of the resulting equation, we arrive at the amplitude-frequency response relation

[TABLE]

There is a sharply increased response when $\Omega$ is close to the resonant frequency associated with the eigenmode, as seen in Figure 5. Different magnitudes of force $F$ are shown. When the force is smaller, the response is much greater, thereby allowing the model to capture a very large range of forcing amplitudes with only relatively small variations in acoustic pressure.

We can also observe (by solving for $\psi_{k}$) that a phase delay of half a cycle is accumulated as we cross the resonant frequency. The group delay, the time required for information to be delivered, is then given by the derivative $-\mathrm{d}\psi/\mathrm{d}\Omega$. It is observed that delays of several milliseconds are observed in the vicinity of resonance (Figure 5).

When studying relations such as (3.24), it becomes apparent that the resonant behaviour occurs slightly away from where is expected (based on the linear system). This is a general property of non-linear systems and can be understood by examining the harmonic response of the unforced non-linear system. Solving (3.24) in the case $F=0$ gives the relationship shown in Figure 6, known as a backbone curve [35]. We see that with increasing amplitude the natural harmonic response of the non-linear system is perturbed away from the resonant frequency.

3.3 Fully-coupled system

Bearing in mind the above analysis of the single-mode approximation (3.7) we now return to the fully-coupled system (3.5). Much of the analysis from Section 3.2 can be readily repeated for the matrix system, particularly with the use of numerical schemes for solving non-linear systems of equations, as was used to produce Figure 3. We focus our attention on the elements which tangibly differ from the above discussions.

If we repeat the eigenvalue analysis of Section 3.2.1 we find that a series of Hopf bifurcations take place, at successive parameter values. When we linearise (3.5) about $\boldsymbol{\alpha}=0$, since the coupling between modes takes place within the non-linear part of the system, we reach $N$ uncoupled linear systems each of which has Jacobian of the form (3.9). This means that each time $\mu$ passes one of the critical values $\mu_{n}^{c}$, $n=1,\dots,N$ (as defined in (3.8)) a Hopf bifurcation occurs.

An important feature of the fully-coupled system, which we highlight since it is not the case for the single-mode formulation, is the ability to predict phase delays of more than half a cycle. In uncoupled oscillator systems the phase delay will not exceed half a cycle (cf. Figure 5) however the cochlea is well-known to exhibit delays of several cycles [62, 3]. Recalling the decomposition (3.2) and the harmonic balance techniques used in Section 3.2.3, the phase delay at a given point $x_{0}$ is given by the complex argument

[TABLE]

An example of how (3.25) varies as a function of the harmonic forcing frequency is shown in Figure 7. The alternating cliffs and plateaus are because the delay increases much more quickly in the region of one of the system’s resonant frequencies.

4 Discussion

We have presented a design for an active acoustic metamaterial that is capable of mimicking the properties of the cochlea. Based on a size-graded array of high-contrast resonators, the structure has similar dimensions to the cochlea and has a resonant spectrum that falls broadly within the range of audible frequencies. This design is able to filter different frequencies in space and, with the introduction of a non-linear amplification term, replicate the fundamental properties of the cochlear amplifier.

A modal decomposition was used to approach the coupled-resonator problem. A thorough analysis of a reduced version of the system was undertaken in order to demonstrate its key properties. The aspects which differ most significantly in the fully-coupled case, such as the accumulation of large phase delays, were accounted for.

A final, interesting, aside is that, while the array of compressible resonators considered in this work is presented as a model for the cochlea on the grounds that it can simulate the oscillations of the basilar membrane, there is evidence that the cochlea itself contains compressible elements [63]. The details of this compressibility, and its relevance to cochlear function, are not yet clear.

The code developed for this study is available online at

https://github.com/davies-b/hopf_active_cochlea.

Acknowledgement

The authors would like to thank Andrew Bell and A. James Hudspeth for insightful comments made on an early version of this manuscript.

Appendix A Multipole expansion method

The layer-potential representation (2.6) reduces the Helmholtz problem (2.4) to finding density functions, $\phi$ and $\psi$, such that the two transmission conditions on $\partial D$ are satisfied. That is, we wish to find $\phi,\psi\in L^{2}(\partial D)$ such that

[TABLE]

where equality holds as elements of $L^{2}(\partial D)$ and the notation $\partial_{\nu}^{-}$ and $\partial_{\nu}^{+}$ denotes the derivative in the direction of the normal to $\partial D$, from the inside and outside of $D$, respectively. More details on the use of layer potentials in solving scattering problems can be found in e.g. [37, 55].

Since we are interested in the case of circular resonators, $\phi$ and $\psi$ are, on each $\partial D_{n}$, $2\pi$-periodic functions of $\theta_{n}$ where $(r_{n},\theta_{n})$ denotes a polar coordinate system about on the centre of $D_{n}$. Such functions admit Fourier expansions of the form

[TABLE]

for coefficients $a_{m}^{n}$, and similarly for $\psi$. The reason such an expansion is useful is that $\mathcal{S}_{D_{n}}^{k}[e^{im\theta_{n}}]$ has an explicit representation, shown in [64] to be given by

[TABLE]

where $J_{m}$ and $H_{m}^{(1)}$ are the Bessel and Hankel functions of the first kind, respectively, $c_{n}=-\frac{i\pi R_{n}}{2}$ and $R_{n}$ is the radius of $D_{n}$.

In order to apply this method to the case of $N\in\mathbb{N}$ resonators we also require an expression for $\mathcal{S}_{D_{n^{\prime}}}^{k}[e^{im\theta_{n}}]$, where $n\neq n^{\prime}$. This is achieved through the use of Graf’s addition formula [55], which says that for any $x,y\in\mathbb{R}^{2}$ such that $|x|>|y|$, it holds that the Helmholtz green’s function $\Gamma^{k}$ is given by

[TABLE]

where $x=(|x|,\vartheta_{x})$ and $y=(|y|,\vartheta_{y})$ are polar representations around a common origin.

Finally, we make the identification $L^{2}(\partial D)\cong L^{2}(\partial D_{1})\times\ldots\times L^{2}(\partial D_{N})$ and decompose the single layer potential as

[TABLE]

where \mathcal{S}_{D_{n^{\prime}}}^{k}\big{|}_{\partial D_{n}}:L^{2}(\partial D_{n^{\prime}})\rightarrow L^{2}(\partial D_{n}) is the evaluation of $\mathcal{S}_{D_{n^{\prime}}}^{k}$ on $\partial D_{n}$. Let $z_{n}^{n^{\prime}}$ be the vector from the centre of $D_{n}$ to that of $D_{n^{\prime}}$, then the off-diagonal terms in (A.4) take the form $\mathcal{S}_{D_{n^{\prime}}}^{k}[e^{im\theta_{n}}]=\mathcal{S}_{z_{n}^{n^{\prime}}+D_{n}}^{k}[e^{im\theta_{n}}]$. The addition of $z_{n}^{n^{\prime}}$ within the integrand can then be decomposed using (A.3).

The derivatives appearing in (A.1) can be handled similarly, based on the expressions

[TABLE]

which can be derived by differentiating (A.2).

Finally, we truncate the Fourier basis on each $\partial D_{n}$, using only $\{e^{im\theta_{n}}:n=-M,\dots,M\}$ for some $M\geq 0$, to reach an approximate matrix representation for (A.1).