2.1. Equations for the bulk

The Navier–Stokes equations for the dimensional fluid velocity $\hat {\boldsymbol{u}}$ and modified pressure $\hat {p}$ in the co-rotating frame are given by (Landau & Lifshitz Reference Landau and Lifshitz1987)

(2.1a) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\hat {\boldsymbol{u}}&=0, \end{align}

(2.1b) \begin{align} \rho \left (\frac {D \boldsymbol{\hat {u}}}{D \hat {t}}+\frac {\partial \boldsymbol{\hat {\varOmega }}}{\partial \hat {t}}\times \boldsymbol{\hat {x}}+2 \boldsymbol{\hat {\varOmega }}\times \boldsymbol{\hat {u}}+\boldsymbol{\hat {\varOmega }}\times (\boldsymbol{\hat {\varOmega }}\times \boldsymbol{\hat {x}})\right )&=-\boldsymbol{\nabla } \hat {p}+\mu {\nabla }^2\boldsymbol{\hat {u}}. \end{align}

Here, the first of the extra terms on the left-hand side is the Euler force, the second term is the Coriolis force and the final additional term is the centrifugal force, each due to the imposed rotation. The pressure $\hat {p}$ is a modified pressure in the sense that it incorporates the linear fictitious force associated with the linear acceleration of the SCC (Buskirk et al. Reference Buskirk, Watts and Liu1976). The fluid is assumed to satisfy the no-slip condition at the edges of the walls, so that $\boldsymbol{\hat {u}}(\hat {r}=\hat {a}(\hat {s}))=\boldsymbol{0}$ for $\hat {s}\in \,(0,2\pi R)$ . Motivated by the small strains in the cupula (Selva et al. Reference Selva, Oman and Stone2009), it is modelled as a linearly elastic material, satisfying the steady Navier equations:

(2.2a) \begin{align}& \rho _s\left (\frac {\partial ^2 \boldsymbol{\hat {u}}_s}{\partial t^2}+\frac {\partial \boldsymbol{\hat {\varOmega }}}{\partial \hat {t}}\times \boldsymbol{\hat {x}}+2\boldsymbol{\hat {\varOmega }}\times \frac {\partial \boldsymbol{\hat {u}}_s}{\partial \hat {t}}+\boldsymbol{\hat {\varOmega }}\times (\boldsymbol{\hat {\varOmega }}\times \boldsymbol{\hat {x}})\right ) =\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat {\tau }}, \end{align}

(2.2b) \begin{align}&\qquad\quad \boldsymbol{\hat {\tau }}=2\mu _s\boldsymbol{\hat {E}}+\lambda _s \textrm{tr}(\boldsymbol{\hat {E}})\boldsymbol{1},\quad 2\boldsymbol{\hat {E}}=\boldsymbol{\nabla }\boldsymbol{\hat {u}}_s+(\boldsymbol{\nabla }\boldsymbol{\hat {u}}_s)^\intercal , \end{align}

alongside a linear Hookean constitutive law relating stress $\boldsymbol{\hat {\tau }}$ and the cupular displacement field $\boldsymbol{\hat {u}}_s(\boldsymbol{x},t)$ . The Lamé constants are $\mu _s=E/(2(1+\nu _s))$ and $\lambda _s=\nu _s E/((1+\nu _s) (1-2\nu _s))$ (Bower Reference Bower2009). Finally, we model the fluid–structure interaction at the cupula–endolymph boundary in the usual way, imposing continuity of velocity and stress (Gkanis & Kumar Reference Gkanis and Kumar2006).

2.2. Numerical simulations

The system of (2.1)–(2.2) was simulated in COMSOL for different values of the Young’s modulus of the cupula and with a Poisson ratio $\nu _s=0.48$ . We impose a simple sinusoidal forcing, given by $\hat { \varOmega }(\hat {t}) = \varOmega _0\sin (2\pi \hat {t}/\mathcal{T})$ , with $\varOmega _0=1\,\textrm{rad/s}$ and $\mathcal{T}=1\textrm{ s}$ . The geometrical parameters are $a = 1.6 \times 10^{-4}$ m, $R=3.2\times 10^{-3}$ m and $t_h = 0.8\times 10^{-4}$ m (Daocai et al. Reference Daocai, Qing, Ximing, Jingzhen, Cheng and Xiangxing2014).

The equations were solved for a solid cupula of finite size, i.e. with no thin cupula assumption. Further details on the numerical simulations are available in Appendix B.1. This includes deformation profiles of the solid cupula (figure 13). Moreover, figure 13 shows that the magnitude of cupula deformation is inversely proportional to Young’s modulus.

Flow profiles on either side of the cupula produced by these numerical simulations are shown in figure 3. Here, a top view of the mid-plane of the flow around the canal is plotted, with colour indicating the speed of the flow, normalised by the maximum speed throughout the flow, and with fast regions coloured red and stagnant regions coloured blue. The cupula appears as the central region and is shown in its deformed configuration. This deflection is imperceptible for all except the $E=10^2\textrm{ Pa}$ case and so the relative magnitude of the cupula deformation is indicated by the grey scale colouring, which shows that it is maximum at the centre. Figure 3 shows the velocity field at time $\hat {t}=0.25$ s, a stage at which transients from the initial condition remain significant. However, it is not necessary to wait for these transients to decay completely to observe the emergence of a distinctly asymmetric velocity profile.

In the panels of figure 3, the Young’s modulus, $E$ , of the cupula increases from left to right. As should be expected, the magnitude of the deflection decreases as the material becomes stiffer. This is confirmed quantitatively in figure 13 of Appendix B.1. Surprisingly, however, we also observe a significant change in flow behaviour: for small values of $E$ , the flow is axially symmetric about the centreline of the tube, but as $E$ increases, the flow transitions, losing its axisymmetry and exhibiting vortical structures. It is worth noting that these plots use values of the Young’s modulus that cover a typical physiological range for soft biological tissues (Goriely Reference Goriely2017). Notwithstanding our remark in the introduction that even lower values have been suggested via indirect methods, a flow transition with physiologically relevant values of $E$ suggests that asymmetric flow may be present in a physiological vestibular system. If so, this would contradict the assumption of axial symmetry that is typical in previous models (Rabbitt & Damiano Reference Rabbitt and Damiano1992; Obrist Reference Obrist2008), and raises interesting questions about what impact such asymmetry might have on the mechanics of balance and rotational sensing. To study this behaviour further, we thus turn now to a theoretical analysis of the governing system.

Figure 3. Cross-section of the velocity fields in the cupula as computed using COMSOL simulations. Results are shown for a range of cupula stiffnesses. Colour represents the relative magnitude of the fluid speed, with red denoting regions in which the flow is fast and blue representing stagnant regions; streamlines are represented by solid black curves. As the stiffness of the cupula increases, a symmetry-breaking of the flow occurs. In particular, for values of the Young’s modulus $E\gt 10^3$ Pa, the flow is usually not axially symmetric. Here, $\hat { \varOmega }(\hat {t}) =\varOmega _0 \sin (2\pi \hat {t}/\mathcal{T})$ and the snapshots are taken at $\hat {t}=0.25$ s, with $\varOmega _0=1$ rad s $^{-1}$ and $\mathcal{T}=1$ s. The geometrical parameters are $a = 1.6 \times 10^{-4}$ m, $R=3.2\times 10^{-3}$ m and $t_h = 0.8\times 10^{-4}$ m. (a) $E = 10^{2}\,\textrm{Pa}$ , (b) $E = 10^{3}\,\textrm{Pa}$ , (c) $E = 10^{4}\,\textrm{Pa}$ and (a) $E = 10^{5}\,\textrm{Pa}$ .

2.3. Theoretical model

To investigate the symmetry breaking observed in the numerical simulations, and to obtain a qualitative understanding of flow and deformation characteristics, in this section, we use asymptotic analysis to derive a reduced order equation for the deformation of the cupula.

To capture the geometry of the semicircular canal, we introduce toroidal coordinates $(\hat {r},\theta , \hat {s})$ , in which $\hat {s}\,\in \,(0,2\pi R)$ is the arc-length along the centreline of the torus. A sketch of this coordinate system is provided in figure 2(b). Cartesian coordinates are related to toroidal coordinates via (Pedley Reference Pedley1980)

(2.3) \begin{equation} \begin{cases} \hat {x} = &(R + \hat {r }\cos {\theta })\cos \left(\dfrac {\hat {s}}{R}\right)\!, \\[8pt]\hat {y} = &(R+\hat {r}\cos \theta )\sin \left(\dfrac {\hat {s}}{R}\right)\!, \\[8pt]\hat {z} = & -\hat {r} \sin \theta . \end{cases} \end{equation}

The negative sign in the last equation ensures that the orthonormal basis vectors $\{\boldsymbol{e}_r, \boldsymbol{e}_\theta , \boldsymbol{e}_s\}$ follow the right-hand rule. We can now rewrite the Navier–Stokes equations (2.1) in component form. Writing the velocity vector as $\boldsymbol{\hat {u}}=\hat {u} \boldsymbol{e}_r +\hat {v}\boldsymbol{e}_\theta +\hat {w}\boldsymbol{e}_s$ , the continuity equation becomes (see Pedley Reference Pedley1980, for example)

(2.4) \begin{equation} \frac {\partial \hat {u}}{\partial \hat {r}}+\frac {\hat {u}}{\hat {r}}+\frac {1}{\hat {r}}\frac {\partial \hat {v}}{\partial \theta }+\frac {1}{h}\frac {\partial \hat {w}}{\partial \hat {s}}-\frac {\hat {v}\sin \theta }{R h}+\frac {\hat {u}\cos \theta }{R h}=0, \end{equation}

where $h=1+\hat {r}\cos (\theta )/R$ is a scale factor. The momentum equations are (see Pedley Reference Pedley1980 for the equations expressed in an inertial frame)

(2.5a) \begin{align} \rho &\left (\frac {\partial \hat {u}}{\partial \hat {t}} +\hat {u}\frac {\partial \hat {u}}{\partial \hat {r}}+\frac {\hat {v}}{\hat {r}}\frac {\partial \hat {u}}{\partial \theta }+\frac {\hat {w}}{h}\frac {\partial \hat {u}}{\partial \hat {s}}-\frac {\hat {v}^2}{\hat {r}}-\frac {\hat {w}^2}{h}\frac {\cos \theta }{R} -2\hat {\varOmega }\hat {w}\cos \theta -\hat {\varOmega }^2Rh\cos \theta \right ) \nonumber \\&= -\frac {\partial \hat {p}}{\partial \hat {r}}+\frac {\mu }{\hat {r}h}\left [\frac {\hat {r}}{h}\frac {\partial }{\partial \hat {s}}\left (\frac {\partial \hat {u}}{\partial \hat {s}}-\frac {\partial }{\partial \hat {r}}(h\hat {w})\right )-\frac {\partial }{\partial \theta }\left (\frac {h}{\hat {r}}\frac {\partial }{\partial \hat {r}}(\hat {r} \hat {v})-\frac {h}{\hat {r}}\frac {\partial \hat {u}}{\partial \theta }\right )\right ]\!, \end{align}

(2.5b) \begin{align} \rho &\left (\frac {\partial \hat {v}}{\partial \hat {t}}+\hat {u}\frac {\partial \hat {v}}{\partial \hat {r}}+\frac {\hat {v}}{\hat {r}}\frac {\partial \hat {v}}{\partial \theta }+\frac {\hat {w}}{h}\frac {\partial \hat {v}}{\partial \hat {s}}+\frac {\hat {u}\hat {v}}{\hat {r}}+\frac {\hat {w}^2}{Rh}\sin \theta +2\hat {\varOmega } \hat {w}\sin \theta +\hat {\varOmega }^2 Rh\sin \theta \right ) \nonumber \\&=-\frac {1}{\hat {r}}\frac {\partial \hat {p}}{\partial \theta }+\frac {\mu }{h}\frac {\partial }{\partial \hat {r}}\left [\frac {h}{\hat {r}}\left (\frac {\partial }{\partial \hat {r}}(\hat {r} \hat {v})-\frac {\partial \hat {u}}{\partial \theta }\right )\right ]-\frac {\mu }{\hat {r}h^2}\left [\frac {\partial }{\partial \theta }\left (h \frac {\partial \hat {w}}{\partial \hat {s}}\right )-\hat {r}\frac {\partial ^2 \hat {v}}{\partial \hat {s}^2}\right ]\!, \end{align}

(2.5c) \begin{align} \rho &\bigg(\frac {\partial \hat {w}}{\partial \hat {t}}+\hat {u}\frac {\partial \hat {w}}{\partial \hat {r}}+\frac {\hat {v}}{\hat {r}}\frac {\partial \hat {w}}{\partial \theta }+\frac {\hat {w}}{h}\frac {\partial \hat {w}}{\partial \hat {s}}+\frac {\hat {u}\hat {w}}{Rh}\cos \theta -\frac {\hat {v}\hat {w}}{Rh}\sin \theta \bigg . \nonumber \\&\bigg .+\frac {d\hat {\varOmega }}{d\hat {t}}(R+\hat {r}\cos \theta )+2\hat {\varOmega } \hat {u}\cos \theta -2\hat {\varOmega } \hat {v}\sin \theta \bigg ) \nonumber \\&= -\frac {1}{h}\frac {\partial \hat {p}}{\partial \hat {s}}+\frac {\mu }{\hat {r}^2}\frac {\partial }{\partial \theta }\left [\frac {1}{h}\frac {\partial }{\partial \theta }(h\hat {w})-\frac {\hat {r}}{h}\frac {\partial \hat {v}}{\partial \hat {s}}\right ] -\frac {\mu }{\hat {r}}\frac {\partial }{\partial \hat {r}}\left [\frac {\hat {r}}{h}\left (\frac {\partial \hat {u}}{\partial \hat {s}}-\frac {\partial }{\partial \hat {r}}(h \hat {w})\right )\right ]\!. \end{align}

We locate the cupula at arc length position $\hat {s}=0$ . As indicated here, our numerical simulations have shown that the deformation of the cupula is small compared with the tube radius, suggesting that strains are small and thus that it is sufficient to use a linear equation. This small strain assumption will be confirmed in § 2.4 through a scaling argument. We write the equations for the solid deformation of the cupula in cylindrical coordinates, as the cupula is thin enough that the curvature of the SCC is not important. The cupula is thus modelled as a solid cylinder $0\leqslant \hat {r}\leqslant \hat {a}(0)$ , $0\leqslant \theta \leqslant 2\pi$ , $-t_h/2\leqslant \hat {z}\leqslant t_h/2$ , with the centre of mass of the cupula ( $\hat {r}=\theta =\hat {z}=0$ ) located at position $\hat {r}=\theta =\hat {s}=0$ (or equivalently, $\hat {s}=2\pi R$ ) in terms of the toroidal coordinates of the fluid problem (see figure 2 a). The thickness of the cupula $t_h$ is not assumed to be small when compared with its radius $\hat {a}(0)$ . The equations for the solid deformation in component form (2.2) are given in Appendix B.2.

We require boundary conditions for both the fluid problem (2.5) and the elastic problem (2.2). The walls of the endolymph are assumed to have a no-slip condition, so that in the rotating frame, $\boldsymbol{\hat {u}}(\hat {r}=\hat {a})=0$ . We require a second set of boundary conditions where the cupula and the endolymph meet. The kinematic boundary condition requires that the cupular velocity and the endolymph velocity at the surface of the cupula must match. As the spatial gradients of the cupula’s deformation are small, we can write this as

(2.6a) \begin{align} &\qquad\quad \frac {\partial \hat {\boldsymbol{u}}_s}{\partial \hat {t}}(\hat {r},\theta ,\hat {z}=t_h/2,\hat t) =\boldsymbol{\hat {u}}(\hat {r}, \hat {s}=t_h/2,\hat {t}), \end{align}

(2.6b) \begin{align} &\frac {\partial \hat {\boldsymbol{u}}_s}{\partial \hat {t}}(\hat {r},\theta ,\hat {z}=2\pi R-t_h/2,\hat t) =\boldsymbol{\hat {u}}(\hat {r},\hat {s}=2\pi R-t_h/2,\hat {t}). \end{align}

We also require boundary conditions for the solid problem (2.2). The precise attachment between the cupula and the utricle is an open area of research and the conditions to be satisfied are not immediately clear. We opt to implement a straightforward choice motivated by Rabbitt & Damiano (Reference Rabbitt and Damiano1992) who model the cupula as a clamped plate, but we extend this to include thickness effects. In this regime, the forcing for the cupular displacement is given by the pressure jump across the two sides of the cupula. In particular, at the flat faces of the cylinder, we impose the following conditions on the stress tensor:

(2.7a) \begin{align}& \hat {\tau }_{zz}=\pm \frac {1}{2}\left [\hat {p}(\hat {r},\theta ,\hat {s}=2\pi R-t_h/2,\hat {t})-\hat {p}(\hat {r},\theta ,\hat {s}=t_h/2,\hat {t})\right ]\equiv \pm \frac {\Delta \hat {p}}{2}, \quad &\text{at }\hat {z}=\pm t_h/2, \end{align}

(2.7b) \begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \hat {\tau }_{rz}=0,\quad \hat {\tau }_{\theta z}=0, \quad &\text{at }\hat {z}=\pm t_h/2. \end{align}

On the curved face, we impose zero displacement, $\boldsymbol{\hat u}_s(\hat {r}=\hat {a}(0))=\boldsymbol 0$ . We also require initial conditions for the velocity, $\boldsymbol{\hat {u}}(\boldsymbol{\hat {x}},\hat {t}=0)$ , and the cupular deflection, $\boldsymbol{\hat {u}_s}(\hat {t}=0)$ . Unless otherwise stated, we assume that the system is initially at rest and is undeformed.

2.4. Scalings and non-dimensionalisation

The SCCs in mammals are thin and slender: with an aspect ratio $\epsilon =a/R$ between 0.05 and 0.1 (Daocai et al. Reference Daocai, Qing, Ximing, Jingzhen, Cheng and Xiangxing2014). It is therefore natural to exploit $\epsilon \ll 1$ and use an asymptotic approach to perform a long wavelength asymptotic analysis similar to lubrication theory. We also know from the continuity equation that if $\hat {w}\sim \mathcal{U}$ , then $\hat {u},\hat {v}\sim \epsilon \mathcal{U}$ . This is required to have a non-trivial balance and is a typical scaling in lubrication theory (Papageorgiou Reference Papageorgiou1995; Craster & Matar Reference Craster and Matar2006). To make this asymptotic intuition more formal, we non-dimensionalise, scaling the dimensional variables according to

(2.8) \begin{align} \begin{split} &\hat {r}=ar,\quad \hat {s} = R s,\quad \hat {z}=az,\quad \hat {t}=\mathcal{T}t,\quad \hat {w} = \frac {R\varOmega _0a^2}{\mathcal{T}\nu }w,\\ & \hat {p} =\frac {\mu R \mathcal{U}}{a^2}p,\quad \hat {\boldsymbol{u}}_s= \frac {a^2 R\varOmega _0}{\nu }\boldsymbol{u}_s,\quad \hat {\boldsymbol{\tau }}=\frac {EaR\varOmega _0}{\nu }\boldsymbol{\tau },\quad \hat {\varOmega }=\varOmega _0\varOmega (t). \end{split} \end{align}

(‘Unhatted’ variables are therefore dimensionless counterparts of the corresponding hatted variables.) Here, $\mathcal{T}$ is the time scale of variation of the forcing and $\nu =\mu /\rho$ is the kinematic viscosity of the endolymph. The velocity scale $\mathcal{U}=R\varOmega _0a^2/(\mathcal T\nu )$ is chosen to balance the viscous forces with the Euler force. The deformation scale is chosen to balance the kinematic condition (2.6a ). The pressure scale is chosen viscously to enforce incompressibility. Note that $\varOmega (t)$ is the dimensionless angular velocity and $\dot {\varOmega }(t)$ is the dimensionless angular acceleration.

2.4.1. Dimensionless equations

Following the rescaling of (2.4)–(2.6a ), the dimensionless continuity equation (2.4) reads

(2.9) \begin{equation} \frac {\partial u}{\partial r}+\frac {u}{r}+\frac {1}{r}\frac {\partial v }{\partial \theta }+\frac {1}{h}\frac {\partial w}{\partial s}-\epsilon \frac {v\sin \theta }{h}+\epsilon \frac {u\cos \theta }{h}=0, \end{equation}

where the scale factor $h$ can be expressed as $h=1 + \epsilon r \cos \theta$ . The dimensionless Navier–Stokes equations (2.5) along the $\boldsymbol{e}_r$ , $\boldsymbol{e}_\theta$ and $\boldsymbol{e}_s$ directions are given respectively by

(2.10a) \begin{align} \epsilon &\textit{St} \frac {\partial u}{\partial t} +\epsilon ^3 \textit{Re}\left ( u\frac {\partial u}{\partial r}+\frac {v}{r}\frac {\partial u}{\partial \theta }+\frac {w}{h}\frac {\partial u}{\partial s}- \frac {v^2}{r}-\frac {1}{\epsilon }\frac {w^2}{h}\cos \theta \right ) \nonumber \\&-2\varOmega _0 \mathcal{T}{St} \varOmega (t)w\cos \theta -\varOmega _0 \mathcal{T} \varOmega (t)^2 h\cos \theta \nonumber \\&= -\frac {1}{\epsilon }\frac {\partial p}{\partial r}+\frac {\epsilon }{rh}\left [\frac {r}{h}\frac {\partial }{\partial s}\left (\epsilon ^2\frac {\partial u}{\partial s}-\frac {\partial }{\partial r}(hw)\right )-\frac {\partial }{\partial \theta }\left (\frac {h}{r}\frac {\partial }{\partial r}(r v)-\frac {h}{r}\frac {\partial u}{\partial \theta }\right )\right ]\!, \end{align}

(2.10b) \begin{align} \epsilon & \textit{St}\frac {\partial v}{\partial t}+\epsilon ^2 \textit{Re}\left (\epsilon u\frac {\partial v}{\partial r}+\epsilon \frac {v}{r}\frac {\partial v}{\partial \theta }+\epsilon \frac {w}{h}\frac {\partial v}{\partial s}+\epsilon \frac {uv}{r}+\frac {w^2}{h}\sin \theta \right ) \nonumber \\&+2\varOmega _0\mathcal{T}{St} \varOmega (t) w\sin \theta +\varOmega _0\mathcal{T}\varOmega (t)^2 h\sin \theta \nonumber \\&=-\frac {1}{\epsilon }\frac {1}{r}\frac {\partial p}{\partial \theta }+\frac {\epsilon }{h}\frac {\partial }{\partial r}\left [\frac {h}{r}\left (\frac {\partial }{\partial r}(r v)-\frac {\partial u}{\partial \theta }\right )\right ]-\frac {\epsilon }{rh^2}\left [\frac {\partial }{\partial \theta }\left (h \frac {\partial w}{\partial s}\right )-\epsilon ^2r\frac {\partial ^2 v}{\partial s^2}\right ]\!, \end{align}

(2.10c) \begin{align} \textit{St}&\frac {\partial w}{\partial t}+\dot {\varOmega }(t)(1+\epsilon r\cos \theta )+\epsilon 2\,\textit{St}\,\mathcal{T}\varOmega _0\varOmega (t)(u\cos \theta -v\sin \theta ) \nonumber \\+&\epsilon ^2 \textit{Re}\left (u\frac {\partial w}{\partial r}+\frac {v}{r}\frac {\partial w}{\partial \theta }+\frac {w}{h}\frac {\partial w}{\partial s}+\epsilon \frac {uw}{h}\cos \theta -\epsilon \frac {vw}{h}\sin \theta \right ) \nonumber \\ &=-\frac {1}{h}\frac {\partial p}{\partial s}+\frac {1}{r}\frac {\partial }{\partial r}\left [\frac {r}{h}\left (\frac {\partial }{\partial r}(h w)-\epsilon ^2\frac {\partial u}{\partial s}\right )\right ]+\frac {1}{r^2}\frac {\partial }{\partial \theta }\left [\frac {1}{h}\frac {\partial }{\partial \theta }(hw)-\epsilon ^2\frac {r}{h}\frac {\partial v}{\partial s}\right]\!. \end{align}

The dimensionless form of (2.2) is given by

(2.11) \begin{align} \epsilon \frac {\rho _s}{\rho }\frac {1}{\kappa }\left ( {St}\frac {\partial ^2\boldsymbol{u}_s}{\partial t^2}+\frac {\partial \boldsymbol{\varOmega }}{\partial t}\times \boldsymbol{x}+2{St}\varOmega _0\mathcal{T}\boldsymbol{\varOmega }\times \frac {\partial \boldsymbol{u}_s}{\partial t}+\varOmega _0\mathcal{T}\boldsymbol{\varOmega }\times (\boldsymbol{\varOmega }\times \boldsymbol{x})\right ) =\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }. \end{align}

The only inhomogeneous boundary conditions for the solid problem are

(2.12) \begin{align} \tau _{zz}(r,\theta ,z=\pm \beta /2,t)=\pm \Delta p(r,\theta ,t)/(2\kappa ), \end{align}

with $\beta =t_h/a$ the dimensionless thickness of the cupula. The non-dimensionalisation procedure introduces several dimensionless parameters. We shall see that the most important of these are the stiffness

(2.13) \begin{equation} \kappa =\frac {E\mathcal{T}a}{R\mu }, \end{equation}

which measures the time scale of forcing to the time scale of the cupula’s relaxation, and the Stokes number

(2.14) \begin{equation} \textit{St}= \frac {a^2}{\nu \mathcal{T}}, \end{equation}

which measures the time scale of vorticity diffusion across the channel width, $a^2/\nu$ , to the time scale of motion, $\mathcal{T}$ . (Note that this version of the Stokes number arises in Stokes’ second problem and is sometimes replaced by the Womersley number, $\textit{Wm}={St}^{1/2}$ .) We also introduce the Reynolds number of the flow, ${\textit{Re}}=\rho \mathcal{U}R/\mu$ , as well as the density ratio, $\rho _s/\rho$ , and the time scale ratio, $\varOmega _0\mathcal{T}$ . To understand the relative size (and hence importance) of these parameters, we next discuss characteristic parameter values and typical sizes of dimensionless parameters.

3.1. Expanded solution

Substitution of the formal expansion equation (3.1) into the governing equation (2.10) yields the following leading order balance:

(3.2a) \begin{align} \dot {\varOmega }(t) &=-\frac {\partial p_0}{\partial s}+\frac {1}{r}\frac {\partial }{\partial r}\left (r\frac {\partial w_0}{\partial r}\right )+\frac {1}{r^2}\frac {\partial ^2w_0}{\partial \theta ^2}, \end{align}

(3.2b) \begin{align} &\qquad\quad\frac {\partial p_0}{\partial r}=\frac {\partial p_0}{\partial \theta }=0, \end{align}

(3.2c) \begin{align} &\frac {\partial u_0}{\partial r}+\frac {u_0}{r}+\frac {1}{r}\frac {\partial v_0}{\partial \theta }+\frac {\partial w_0}{\partial s}=0, \end{align}

(3.2d) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }_0=&\boldsymbol 0,\quad \tau _{zz0}(r,\theta ,z=\pm \beta /2,t)=\pm \frac {1}{2\kappa }\Delta p_0(t). \end{align}

We identify the usual lubrication/boundary layer theory result that the leading order pressure is constant along a cross-section (Craster & Matar Reference Craster and Matar2006; Tavakol et al. Reference Tavakol, Froehlicher, Holmes and Stone2017). This has the direct consequence that the pressure difference $\Delta p_0(t)=\int _{\epsilon \beta /2}^{2\pi -\epsilon \beta /2}\partial p_0/\partial s\,\textrm {d} s$ is only a function of time. Hence, we may seek a solution to the solid problem of the form $\boldsymbol{u}_{s0}(r,\theta ,z,t)=\Delta p_0(t)\boldsymbol{\bar {u}}_{s0}(r,z)/\kappa$ , where $\boldsymbol{\bar {u}}_{s0}=(\bar {u}_{s0}(r,z),0,\bar {w}_{s0}(r,z))^\intercal$ is an axisymmetric solution to (2.11) satisfying the normalised boundary condition $\bar {\tau }_{zz0}=\pm 1/2$ at $z=\pm \beta /2$ . We can thus solve the leading order solid problem independently from the fluid problem. Moreover, motivated by results from the numerical simulations indicating $w_{s0}(r,z)$ varies slowly over $z$ , we introduce the depth-averaged cupular deflection $\eta (r,\theta ,t)=\eta _0(r,t)+\epsilon \,\eta _1(r,\theta ,t)+{\cdots}$ defined as

(3.3) \begin{align} \eta (r,\theta ,t)=\frac {1}{\beta }\int _{-\beta /2}^{\beta /2}w_s(r,\theta ,z,t)\,\textrm {d} z. \end{align}

In Appendix B.2, we show that when ${St}\ll 1$ , we can decompose $\eta$ multiplicatively as

(3.4) \begin{align} \eta _0(r,\theta ,t)=\frac {\Delta p_0(t)}{\kappa }\bar {\eta }(r,\theta ),\quad \bar {\eta }(r,\theta )=\frac {1}{\beta }\int _{-\beta /2}^{\beta /2}\bar {w}_s(r,z,t)\,\textrm {d} z, \end{align}

so that, to leading order, the cupular deflection and the pressure jump are proportional. We may obtain a polynomial approximation for $\bar {\eta }_0(r)$ using techniques from Barber (Reference Barber2010). The final solution (see Appendix B.2 for details on the calculation) is

(3.5) \begin{align} \bar {\eta }_0(r;\beta )=\frac {1}{\beta }\!\int _{-\beta /2}^{\beta /2}\! w_{s0}(r,z;\beta )\,\textrm {d} z=\frac {3}{16}\frac {1-\nu _s^2}{\beta ^3}(1-r^2)^2+\frac {1}{20}\frac {(1+\nu _s)(12-\nu _s)}{\beta }(1-r^2). \end{align}

Having obtained a form of solution for the solid problem in terms of the pressure jump, we turn our attention to the fluid problem. We remark that the differential operators acting on the continuity and momentum equations are the same as in cylindrical coordinates (Batchelor Reference Batchelor1973), and invoking symmetry, we now seek an axisymmetric solution with $v_0=0$ and leading order terms independent of $\theta$ , with $w_0=w_0(r,s,t)$ . The leading order kinematic boundary condition may be written as

(3.6) \begin{align} \frac {\partial \eta _0}{\partial t}=w_0(r,s=0,t)=w_0(r,s=2\pi ,t),\quad u_0(r,s=0,t)=u_0(r,s=2\pi ,t)=0. \end{align}

Our starting point is the continuity equation, which can be integrated over a cross-section to deduce that the flux $Q = \int _0^{2\pi }\int _0^{a(s)}rw(r,\theta , s, t)\,\textrm {d} r\,\textrm {d} \theta$ is conserved in the $s$ direction, i.e. $\partial Q/\partial s=0$ . This means that the flux is exclusively a function of time, a fact we will exploit to derive a reduced-order equation. Turning our attention to the $\mathcal{O}(\epsilon )$ problem, the equations read

(3.7a) \begin{align} \dot {\varOmega }(t) r\cos \theta &=-\frac {\partial p_1}{\partial s}+ r\cos \theta \frac {\partial p_0}{\partial s}+\frac {1}{r}\frac {\partial }{\partial r}\left (r\frac {\partial w_1}{\partial r}\right )+\frac {1}{r^2}\frac {\partial ^2w_1}{\partial \theta ^2}+\cos \theta \frac {\partial w_0}{\partial r}, \end{align}

(3.7b) \begin{align} &\frac {\partial p_1}{\partial r}=\varOmega _0\mathcal{T}\varOmega (t)^2\cos \theta ,\quad \frac {1}{r}\frac {\partial p_1}{\partial \theta }=-\varOmega _0\mathcal{T}\varOmega (t)^2\sin \theta , \end{align}

(3.7c) \begin{align} &\quad\frac {\partial u_1}{\partial r}+\frac {u_1}{r}+\frac {1}{r}\frac {\partial v_1}{\partial \theta }+\frac {\partial w_1}{\partial s}=\cos \theta \left (r\frac {\partial w_0}{\partial s}-u_0\right )\!, \end{align}

(3.7d) \begin{align} &\quad\,\,\frac {\rho _s}{\rho }\frac {1}{\kappa }\left (\frac {\partial \boldsymbol{\varOmega }}{\partial t}\times \boldsymbol{x} +\varOmega _0\mathcal{T}\boldsymbol{\varOmega }\times (\boldsymbol{\varOmega }\times \boldsymbol{x})\right )=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }_1, \end{align}

(3.7e) \begin{align} &\qquad\quad\tau _{zz1}(r,z=\pm \beta /2,t)=\pm \Delta p_1(r,t)/(2\kappa ), \end{align}

(3.7f) \begin{align} &\qquad\qquad \Delta p_1(r,t)=\Delta p_1^{\textit{outer}}(t)+\Delta p_1^{\textit{BL}}(r,t). \end{align}

Here, the pressure difference due to the outer flow is $\Delta p_1^{\textit{outer}}(t)=\int _{\epsilon \beta /2}^{2\pi -\epsilon \beta /2}\partial p_1/\partial s\,\textrm {d} s$ and $\Delta p_1^{\textit{BL}}(r,t)$ is a contribution from the boundary layer that forms near the cupula (see Appendix F for details). We note that $p_1(r,\theta , s,t)$ may be decomposed into a pressure gradient along the duct axis and an $s$ -independent pressure variation due to centrifugal effects, so that

(3.8) \begin{equation} p_1(r,\theta ,s,t) = \bar {p}_1(s,t)+\varOmega _0\mathcal{T}\varOmega (t)^2r\cos \theta , \end{equation}

and as we only require the $s$ component of the pressure gradient in the computation of the axial velocity, we may safely ignore the $s$ -independent component of $p_1$ and use $\bar {p}_1(s,t)$ in its place. Moreover, the first-order pressure jump across the cupula is $\Delta p_1^{\textit{outer}}=p_1(s=2\pi ,t)-p_1(s=0,t)=\Delta \bar {p}_1^{\textit{outer}}$ .

The solutions giving the first two orders of the axial velocity in terms of the pressure gradients may be determined directly: $w_0$ is found from (3.2a ), under the assumption of axisymmetry, while $w_1$ may be found by decomposing it into axisymmetric and asymmetric parts, and using standard methods. We find that

(3.9a) \begin{align} &\qquad\qquad\quad w_0(r,s,t)=-\frac {1}{4}\left (\dot {\varOmega }(t)+\frac {\partial p_0}{\partial s}\right )(a(s)^2-r^2), \end{align}

(3.9b) \begin{align} &w_1(r,\theta ,s,t)= -\frac {1}{4}\frac {\partial \bar {p}_1}{\partial s}(a(s)^2-r^2)-\frac {1}{16}\left (\dot {\varOmega }(t)-3\frac {\partial p_0}{\partial s}\right )r(a(s)^2-r^2)\cos \theta . \end{align}

Figure 4. Velocity profiles predicted by (3.9) as the torus aspect ratio, $\epsilon$ , and leading-order term vary. Note how when $\dot {\varOmega }(t)$ and $\kappa \Delta p/(2\pi )$ cancel each other, the asymmetric flow dominates. Moreover, the symmetry breaking becomes observable earlier for larger values of $\epsilon$ , though we reiterate that our theory is formally valid only for $\epsilon \ll 1$ . Note that the horizontal axis may be interpreted as the phase difference between $\dot \varOmega$ and $-\Delta p$ .

We can get a first hint of the symmetry breaking mechanism observed in figure 3 by considering when the asymmetric correction term $\epsilon w_1$ is of a similar size as the symmetric leading order solution $w_0$ . Indeed, it is easy to verify that this will be the case when the pressure gradient approximately cancels out the forcing $\dot \varOmega (t)$ , such that the modulating coefficient $\dot {\varOmega }(t) + ({\partial p_0}/{\partial s})$ in (3.9a ) is close to zero. This is visualised in figure 4, where we plot the velocity $w=w_0+\epsilon w_1$ for several values of $\epsilon$ and $\dot {\varOmega }(t)+ ({\partial p_0}/{\partial s})$ , observing an asymmetric profile when $\dot {\varOmega }(t) + ({\partial p_0}/{\partial s})\ll 1$ . We note that the above-mentioned solutions are only valid far from the cupula, in the slender regions of the SCC ( $s\gg \epsilon )$ , and thus we refer to (3.9) as the outer solution. Closer to the cupula, a boundary layer develops due to the leading order velocity having to adjust from the outer solution $w_0(r,s,t)$ to the cupular deflection profile $\partial \eta _0/\partial t$ as imposed by (3.6) (see the curved streamlines close to the cupula in figure 3). The velocity adjustment causes the aforementioned pressure jump $\Delta p_1^{\textit{BL}}(r,t)$ , which we compute in Appendix F using a matched asymptotic expansion:

(3.10) \begin{align} \begin{split} \Delta p_1^{\textit{BL}}(r,t)&=-\frac {\pi }{I_4}\left (\dot {\varOmega }(t)+\frac {\Delta p_0(t)}{2\pi }\right )f_3(\beta ){\textrm{Re}}\left \{\sum _{n=1}^\infty \bar {A}_nJ_0(\mu _nr)\right \},\\ f_3(\beta )&=\frac {5(1-\nu _s)}{4\beta ^2(12-\nu _s)+10(1-\nu _s)}, \end{split} \end{align}

where $\bar {A}_n\in \mathbb{C}$ are constants independent of time and cupular thickness, and $\mu _n\in \mathbb{C}$ are the solutions in the first quadrant of $J_1(z)^2=J_0(z)J_2(z)$ (for details, see Davis Reference Davis1990). From the velocities (3.9), we may compute the flux, noting that the asymmetric components of $w_1(r,\theta ,s,t)$ integrate to zero because of the $\cos \theta$ term:

(3.11a) \begin{align}& Q_0 = -\frac {2\pi }{16}\left (\dot {\varOmega }(t)+\frac {\partial p_0}{\partial s}\right )a(s)^4, \end{align}

(3.11b) \begin{align}&\qquad Q_1 = -\frac {2\pi }{16}\frac {\partial \bar {p}_1}{\partial s}a(s)^4. \end{align}

Since the flux $Q=Q_0+\epsilon Q_1+\ldots$ is independent of $s$ , we can now integrate the above equations along the axis of the duct, obtaining

(3.12a) \begin{align}& I_4 Q_0 =-\frac {\pi }{8}\left (2\pi \dot {\varOmega }(t)+\Delta p_0\right )\!, \end{align}

(3.12b) \begin{align}&\qquad I_4 Q_1 = -\frac {\pi }{8}\Delta p_1^{\textit{outer}}, \end{align}

where we have defined $I_4=\int _0^{2\pi }a(s)^{-4}\,\textrm {d} s$ and used the fact that $\Delta p_1^{\textit{outer}}=\Delta \bar {p}_1^{\textit{outer}}$ . We remark that when the thickness of the cupula is appreciable, the integrals should be performed from $s=\epsilon \beta /2$ to $s=2\pi -\epsilon \beta /2$ , leading to a slight modification of (3.12), as discussed in Appendix E.

To connect the flow to the cupula displacement, we evaluate the flux using the velocity $w$ at the cupula, where the fluid velocity must equal the velocity of the cupula: $w = ({\partial \eta }/{\partial t})$ . Therefore, we may write $Q_i=\int _0^{2\pi } \int _0^{a_0}r ({\partial \eta _i}/{\partial t})\,\textrm {d} r\,\textrm {d} \theta$ , where $a_0=a(0)$ . In Appendix B.2, we show that the first-order solid problem may be solved as

(3.13) \begin{align} \eta _1(r,\theta ,t)& =\frac {\Delta p_1^{\textit{outer}}(t)}{\kappa }\bar {\eta }_0(r)+\frac {\rho _s}{\rho \kappa }\dot {\varOmega }(t)\bar {\eta }^{\textit{Euler}}_1(r) \nonumber\\ & \quad +\frac {\rho _s\varOmega _0\mathcal{T}}{\rho \kappa }\varOmega (t)^2\bar {\eta }^{\textit{centrif}}_1(r)\cos \theta -\frac {\pi }{I_4\kappa }\left (\dot {\varOmega }(t)+\frac {\Delta p_0(t)}{2\pi }\right )f_3(\beta )\bar {\eta }^{\textit{BL}}_1(r), \end{align}

where $\bar {\eta }_0$ was defined in (3.5) and other $\bar {\eta }_1$ are given in Appendix B.2. Computing the flux, we find

(3.14a) \begin{align} &Q_0=\frac {1}{\kappa }\frac {\textrm {d} \Delta p_0}{\textrm {d} t}\int _0^{2\pi }\int _0^{a(0)}r\bar {\eta }_0(r;\beta )\,\textrm {d} r\,\textrm {d} \theta =\frac {2\pi \alpha _0(\beta )}{\kappa }\frac {\textrm {d} \Delta p_0}{\textrm {d} t}, \end{align}

(3.14b) \begin{align} &Q_1= \frac {2\pi }{\kappa }\left [ \frac {\textrm {d} \Delta p_1^{\textit{outer}}}{\textrm {d} t}\alpha _0(\beta )+\frac {\rho _s}{\rho }\ddot {\varOmega }(t)\alpha _1^{\textit{euler}}(\beta )-\frac {\pi }{I_4}\left (\ddot {\varOmega }(t)+\frac {1}{2\pi }\frac {\textrm {d} \Delta p_0}{\textrm {d} t}\right )f_3(\beta )\alpha _1^{\textit{BL}}(\beta )\!\right ]\!, \end{align}

where the $\alpha$ factors are

(3.15) \begin{align} \alpha _0(\beta )\equiv \int _0^{a(0)}r\bar {\eta }_0(r;\beta )\,\textrm {d} r=\frac {1-\nu _s^2}{32\beta ^3}+\frac {(12-\nu _s)(1+\nu _s )}{80\beta }, \end{align}

and similarly for the other contributions, with explicit forms given in Appendix B.2. The centrifugal term integrates to zero, and we substitute the flux terms $Q_0$ , $Q_1$ into (3.12) from which we obtain a pair of ordinary differential equations for $\Delta p(t)$ :

(3.16a) \begin{align} \frac {\alpha _0(\beta )}{\kappa }\frac {\textrm {d} \Delta p_0}{\textrm {d} t}&=-\frac {1}{16 I_4}\left (2\pi \dot {\varOmega }(t)+ \Delta p_0\right )\!, \end{align}

(3.16b) \begin{align} \frac {\alpha _0(\beta )}{\kappa }\frac {\textrm {d} \Delta p_1^{\textit{outer}}}{\textrm {d} t}&=-\frac {1}{16 I_4} \Delta p_1^{\textit{outer}}-\frac {\rho _s}{\rho }\frac {\alpha _1^{\textit{Euler}}(\beta )}{\kappa }\ddot {\varOmega }(t)\nonumber\\&\quad +\frac {\pi }{I_4\kappa }\alpha _1^{\textit{BL}}(\beta )f_3(\beta )\frac {\textrm {d} }{\textrm {d} t}\left (\dot {\varOmega }(t)+\frac {\Delta p_0}{2\pi }\right ). \end{align}

The boundary layer forcing term can be written in terms of $\Delta \ddot {p}_0$ using (3.16a ), but it will be more useful to leave it in terms of $\Delta p_0$ . We may solve (3.16a ) for any forcing via the integral

(3.17) \begin{equation} \Delta p_0 = -\frac {\pi \kappa }{8\alpha _0(\beta )I_4}\int _0^t\dot {\varOmega }(\tau )e^{-\frac {\kappa }{16I_4 \alpha _0(\beta )}(t-\tau )}\,\textrm {d} \tau . \end{equation}

3.1.1. Expressions for the velocities

Once the pressure jump is known from (3.16a ), the pressure gradient may be computed from

(3.18) \begin{equation} \frac {\partial p_0}{\partial s}=\frac {\Delta p_0}{I_4a(s)^4}+\dot {\varOmega }(t)\left [\frac {2\pi }{I_4a(s)^4}-1\right ],\quad \frac {\partial \bar {p}_1}{\partial s}=\frac {\Delta p_1^{\textit{outer}}}{I_4 a(s)^4}. \end{equation}

Substitution into the axial velocity (3.9) then yields

(3.19a) \begin{align} &w_0 =\frac {\pi }{2I_4a(s)^4}\left [\varOmega (t)+\frac {\Delta p_0}{2\pi }\right ]\left[r^2-a(s)^2\right]\!, \end{align}

(3.19b) \begin{align} &w_1 = \frac {r\left(r^2-a(s)^2\right)\cos \theta }{16I_4 a(s)^4}\left [-3\Delta p_0 +\big(-6\pi +4I_4 a(s)^4\big)\varOmega (t)\!\right ]+\frac {\Delta p_1^{\textit{outer}}}{4I_4a(s)^4} \big(r^2-a(s)^2 \big). \end{align}

The radial and azimuthal velocities can be recovered from the continuity equation,

(3.20a) \begin{align} u_0 &= \frac {\pi }{2I_4a(s)^5}\left [\varOmega (t)+\frac {\Delta p_0}{2\pi }\right ]r\big[r^2-a(s)^2\big]\frac {\textrm {d} a}{\textrm {d} s}, \end{align}

(3.20b) \begin{align} v_0&=0, \end{align}

(3.20c) \begin{align} u_1 & = \frac { 2 \varOmega (t) \left [2 I_4 a(s)^4+\pi \right ]+ \Delta p_0}{16 I_4 a(s)^5}r^2 \big [r^2-a(s)^2\big ] \cos \theta \frac {\textrm {d} a}{\textrm {d} s}+\frac {\Delta p_1^{\textit{outer}}}{4I_4a(s)^5}r(r^2-a(s)^2)\frac {\textrm {d} a}{\textrm {d} s}, \end{align}

(3.20d) \begin{align} v_1 & = -\frac {5 \left [2 \varOmega (t) \left (2 I_4 a(s)^4+\pi \right )+ \Delta p_0\right ]}{16 I_4 a(s)^5}r^2\big [r^2-a(s)^2\big ] \sin \theta \frac {\textrm {d} a}{\textrm {d} s}. \\[9pt] \nonumber \end{align}

Note that for a perfect, i.e. uniform, torus, $a'(s)=0$ and there are no radial or azimuthal velocities.

3.2. Deformation regimes

For a fixed aspect ratio, the key remaining dimensionless parameter is the relative stiffness $\kappa$ (which can be varied by changing the bending stiffness or time scale of forcing). We shall see that changing $\kappa$ makes the system sensitive either to the angular acceleration or the angular velocity. In this section, we analyse how these limiting cases arise by solving (3.16a ). Although (3.16a ) can be solved analytically for any forcing through (3.17), qualitative information may be obtained by considering the large and small $\kappa$ limits.

When the cupula is soft ( $\kappa \ll 1$ ), the solution to (3.16a ) is approximately given by

(3.21) \begin{align} \Delta p_0 \sim -\frac {\pi \kappa }{8\alpha _0(\beta ) I_4}\int _0^t \dot {\varOmega }(\tau )\,\textrm {d} \tau = -\frac {\pi \kappa }{8I_4\alpha _0(\beta )}\varOmega (t) \quad \text{as}\, \kappa \to 0. \end{align}

At the other extreme, for a stiff cupula, characterised by $\kappa \gg 1$ , an approximate solution of (3.17a ) is given by

(3.22) \begin{align} \Delta p_0 \sim -2\pi \dot {\varOmega }(t)+\frac {32I_4\pi \alpha _0(\beta )}{\kappa }\ddot {\varOmega }\quad \text{as}\, \kappa \to \infty . \end{align}

In physical terms, the two results in (3.21) and (3.22) represent qualitatively distinct regimes for the response of the cupular displacement to the imposed rotation of the canal: for a soft cupula ( $\kappa \ll 1$ ), the pressure difference across the cupula, and thus the cupula deformation, is proportional to the angular velocity of the imposed rotation, $\varOmega (t)$ , while for a stiff cupula ( $\kappa \gg 1$ ), the deformation instead follows the angular acceleration $\dot {\varOmega (t)}$ . Since the cupular deformation is thought to be what is detected by the nerve cells in the cupula, this suggests that the cupula can detect either the angular velocity or the angular acceleration to which it is subject – depending on the value of $\kappa$ .

In the latter case, it is easy to verify that the leading order radial and axial velocities suffer a cancellation, as the prefactor of the leading order velocity (3.9a ) is $\dot {\varOmega }(t) + ({\Delta p_0}/{2\pi}) = \mathcal{O} ({1}/{\kappa } )$ , and the asymmetric order $\epsilon$ correction $w_1(r,\theta ,s,t)$ dominates for any $\epsilon$ provided $\kappa$ is sufficiently large (the symmetric component of $w_1$ also scales as $\mathcal{O}(\kappa ^{-1}))$ . This cancellation accounts for the behaviour observed in figure 3: as the Young’s modulus is increased, making the cupula stiffer to the point that $\kappa \gg 1$ , the flow ceases to be symmetric. This is intuitive, as a completely rigid cupula does not allow a net flux and hence the axisymmetric leading order flow must vanish.

We emphasise that although symmetry breaking arises from the breakdown of the asymptotic ordering between the first and second terms in the series, this does not imply a loss of asymptotic ordering in the higher-order terms. The symmetry breaking results from a catastrophic cancellation in the leading-order term, while the first correction remains $\mathcal{O}(\epsilon )$ , and the subsequent terms are expected to retain their anticipated scaling – indicating that the series remains well behaved. We confirm this in the next section by numerically solving the full nonlinear problem (2.1).

To estimate the critical value of $\kappa$ where the transition occurs, we may write $\Delta p_0 =Ae^{2\pi i t} + \textrm{c.c.}$ and $\dot {\varOmega }(t)=B e^{2\pi it} + \textrm{c.c.}$ , with $A,B\in \mathbb{C}$ , and seek the range of $\kappa$ for which $\Delta p$ is in phase with $\varOmega$ . Direct substitution into (3.16a ) yields

(3.23) \begin{equation} A=-\frac {2\pi B}{32 I_4\pi \alpha (\beta )i/\kappa +1}. \end{equation}

Therefore, the transition change occurs at $\kappa _c = 32 \pi I_4 \alpha (b)$ . For $\kappa \lt \kappa _c$ , the phase difference between $A$ and $B$ is more than $({\pi}/{4})$ , and it will be less than $ ({\pi}/{4})$ for $\kappa \gt \kappa _c$ . When the phase difference is small, the response is roughly in phase with the forcing $\varOmega (t)$ , and vice versa.

3.2.1. Sample head rotation

To visualise the different regimes we have identified above, we solve the equation for the pressure jump (3.16a ) with a specific choice for the forcing $\dot \varOmega (t)$ . While many forms of $\dot \varOmega (t)$ could be considered, we are motivated by a clinical head manoeuvre that may be thought of as modelling a slow rotation of the head from right to left. Although other clinical models are described by high-order polynomials (Boselli et al. Reference Boselli, Obrist and Kleiser2009, Reference Boselli, Obrist and Kleiser2013), we choose a particular form that facilitates an analytical solution of the equations, namely

(3.24) \begin{equation} \dot {\varOmega }(t)=\begin{cases} \sin 2\pi t & t\in (0,1),\\ 0 & t\gt 1. \end{cases} \end{equation}

Solving for the pressure gradient through (3.17) yields

(3.25) \begin{align} -\Delta p_0 &= 2\pi \gamma \kappa \int _0^t\varOmega (\tau )e^{-\gamma \kappa (t-\tau )}\,\textrm {d} \tau \nonumber\\&=\begin{cases}2\pi \gamma \kappa \dfrac {2 \pi e^{-\gamma \kappa t}+\gamma \kappa \sin 2 \pi t-2 \pi \cos 2 \pi t}{\gamma ^2 \kappa ^2+4 \pi ^2}& t\lt 1,\\[8pt] \dfrac {4\pi ^2\gamma \kappa }{4\pi ^2+\gamma ^2\kappa ^2}e^{-\gamma \kappa t}(1-e^{-\gamma \kappa })&t\gt 1. \end{cases} \end{align}

Here, $\gamma = 1/[16\alpha _0(\beta )I_4]$ accounts for domain irregularities and cupular thickness. The cupular deformation, which in this regime is proportional to the pressure jump, is plotted in figure 5(a) for different values of $\kappa$ , from which we can clearly see the transition from $\Delta p$ tracking the angular velocity $\varOmega (t)$ for small $\kappa$ to tracking the angular acceleration $\dot {\varOmega }(t)$ for large $\kappa$ , as expected from the preceding analysis. There is an interesting transition region for $\kappa \sim 1$ , where we can see an ‘overshoot’ region at the end of the manoeuvre that has not decayed. This is not the case for either of the limiting regions, where the pressure jump (and cupular displacement) is identically zero after the completion of the head turn. In figure 5(b), we compare how similar the response is to either the angular velocity or the angular acceleration by computing the correlation between the respective functions; for two functions $f(t)$ and $g(t)$ , this correlation is defined as

(3.26) \begin{equation} R(f,g)=\frac {\int _0^T f(t)g(t)\,\textrm {d} t}{\left (\int _0^T f(t)^2\,\textrm {d} t\boldsymbol{\cdot }\int _0^T g(t)^2\,\textrm {d} t\right )^{1/2}}. \end{equation}

As expected from our asymptotic analysis, $\Delta p$ correlates with the angular velocity for small to moderate values of $\kappa$ and the angular acceleration for large values of $\kappa$ . For the parameters used in figure 5, we compute $\kappa _c=32\pi \approx 100$ , in agreement with the transition point observed in the plot.

Figure 5. Influence of dimensionless stiffness $\kappa$ on the cupular deformation. (a) As $\kappa$ is increased, the deformation (normalised by the maximum) transitions from following the angular velocity to following the angular acceleration. (b) This transition with $\kappa$ may be shown by plotting the correlation, $R$ , between the deformation and the angular velocity $\varOmega (t)$ (solid curve) or the angular acceleration $\dot {\varOmega }(t)$ (dashed curve). A transition between the two regimes occurs at $\kappa \approx 100$ . In both plots, colour is used to show the value of $\kappa$ , as indicated in the inset of (a).

5.1. A simple example

For the special case when $a(s)\equiv 1$ , i.e. the tube radius is constant, the pressure gradient can be assumed to be independent of $s$ , that is, $ ({\partial p_0}/{\partial s}) = ({\Delta p_0}/{2\pi})$ , and (5.4a ) simplifies to

(5.6) \begin{equation} \frac {\alpha _0(\beta )}{\kappa }\frac {\textrm {d} \Delta p_0}{\textrm {d} t}=-2 \sum _{n=1}^\infty \left [\lambda _n^{-2}\int _0^t\left (\dot {\varOmega }(\tau )+ \frac {\Delta p_0}{2\pi }\right )\mathcal{K}_n(t-\tau ;{St})\,\textrm {d}\tau \right ]. \end{equation}

To transform this equation into a more manageable form, we define a complete kernel $\mathcal{K}(x;{St})={St}^{-1}\sum _{n=1}^\infty \lambda _n^{-2}e^{-\lambda _n^2x/{St}}$ .

(5.7) \begin{align} \frac {\alpha _0(\beta )}{\kappa }\frac {\textrm {d} \Delta p_0}{\textrm {d} t}=-2\bar {\dot {\varOmega }}(t)-\frac {1}{\pi }\int _0^t\Delta p_0(\tau )\mathcal{K}(t-\tau ;{St})\,\textrm {d} \tau , \end{align}

where we have introduced $\bar {\dot {\varOmega }}(t)=\int _0^t\dot {\varOmega }(\tau )\mathcal{K}(t-\tau ;{St})\,\textrm {d} \tau$ . Equation (5.7) may be efficiently solved numerically by truncating the infinite series in the kernel and transforming the integral equation into a system of ordinary differential equations (ODEs). This is a standard calculation, with details given in Appendix C.

Figure 8. (a) Solution to (5.7), when $\dot \varOmega (t)=0$ and $\Delta p_0 (t=0)=1$ for different values of the Stokes number, showing underdamped dynamics for large enough ${St}$ . (b) Bifurcation diagram, showing the evolution of ${\textrm{Re}}(\bar {\omega })$ (blue) and $\textrm{Im}(\bar {\omega })$ (red). Markers represent the numerically obtained solution from (5.8) and dashed lines the analytical approximation (5.9). (c) Bifurcation diagram for $\dot \varOmega (t)\sim e^{\textit{it}}$ , showing how the transition between the three regimes depends on both $\kappa$ and ${St}$ . Colour represents the complex angle of $\chi$ , with purple representing $\arg (\chi )=0$ , green is $\arg (\chi )=\pi /2$ and yellow is $\arg (\chi )=\pi$ , as described in the colourbar to the right.

5.2. Fluid inertia can make the cupula underdamped

To understand the effect of inertia in the cupular response, we first consider the case of a cupula that is initially stretched by a pressure jump $\Delta p_0(t=0)=1$ in a frame rotating at constant speed, so that $\varOmega (t)=0$ .

The numerical solution to (5.7) for a range of Stokes numbers is given in figure 8(a). Notice that for sufficiently large ${St}$ , $\Delta p(t)$ exhibits decaying oscillatory behaviour, meaning that the cupula is underdamped; this is in contrast to smaller values of ${St}$ , in which the cupula dynamics show an exponential decay whose rate of decay increases with ${St}$ . To understand this change in behaviour as ${St}$ is increased, we seek an exponential ansatz to solve (5.7), with $\Delta p_0 \sim e^{-\omega t}$ in the simple case when $\dot \varOmega (t)=0$ . Direct substitution yields the following condition:

(5.8) \begin{equation} \frac {\alpha _0(\beta ) \pi }{\kappa {St}}\bar {\omega }=\sum _{n=1}^\infty \frac {1}{\lambda _n^2}\frac {1}{\lambda _n^2-\bar {\omega }}, \end{equation}

where the rescaled growth rate is $\bar {\omega }={St}\,\omega$ . Equation (5.8) provides an equation for $\bar {\omega }$ depending on the parameter $q=\alpha _0(\beta )\pi /(\kappa {St})=\alpha _0(\beta ) \pi \rho \nu ^2R/(Ea^3)$ , which is independent of the forcing time scale $\mathcal{T}$ . Analytical progress can be made by truncating the sum at $n=1$ , i.e. considering only the first term. This leads to

(5.9) \begin{equation} \bar {\omega }=\frac {\lambda _1^2}{2}\left [1\pm \left (1-\frac {4}{q\lambda _1^6}\right )^{1/2}\right ]. \end{equation}

From this, we can identify the critical value for the parameter $q$ at which the transition to underdamped dynamics occurs: $q_c = 4/(\lambda _1)^6\approx 0.0207$ .

A natural question to ask now is if fluid inertia can alter the development of the two flow regimes outlined previously in § 3.2? Proceeding as before, we assume a forcing $\dot {\varOmega }(t)=B e^{ i t}$ and try an ansatz $\Delta p_0 = Ae^{ it}$ . Substitution into the integral (5.7) and neglecting contributions from the initial conditions yields

(5.10) \begin{equation} \frac {i\alpha _0(\beta ) A}{\kappa } = -\sum _{n=1}^\infty \frac {2B+A/\pi }{\lambda _n^2(i{St}+\lambda _n^2)}. \end{equation}

It is convenient to define the function $\mathcal{G}:\mathbb{R}\rightarrow \mathbb{C}$ , given by $\mathcal{G}({St})=\sum _{n=1}^\infty \lambda _n^{-2}/( i{St}+\lambda _n^2)$ . We find that the response (characterised by $A$ ) is related to the forcing (characterised by $B$ ) through

(5.11) \begin{equation} A=\frac {-2B \mathcal{G}({St})}{ i\alpha _0(\beta )/\kappa +\mathcal{G}({St})/ \pi }. \end{equation}

Therefore, we see that the angle of the complex quantity

(5.12) \begin{equation} \chi =\frac {-\mathcal{G}({St})}{\pi i\alpha _0(\beta ) +\kappa \mathcal{G}({St})}=\frac {-1}{\kappa +\pi i\alpha _0(\beta )/\mathcal{G}({St})} \end{equation}

will determine if the deformation follows the angular velocity (if the angle is close to $\pi /2$ ) or the angular acceleration (when the angle is close to $0$ or multiples of $\pi$ ). To achieve analytical progress, we truncate the sum in $\mathcal{G}$ , keeping only the first term, and we find

(5.13) \begin{eqnarray} \frac {{\textrm{Re}}(\chi )}{\textrm{Im}(\chi )}=-\frac {\kappa -\pi \alpha _0(\beta )\lambda _1^2{St}}{\pi \lambda _1^4\alpha _0(\beta )}. \end{eqnarray}

Therefore, the curve in the parameter space $(\kappa , {St})$ separating the two regimes is given implicitly by

(5.14) \begin{eqnarray} \left \vert \frac {\kappa }{\pi \lambda _1^4\alpha _0(\beta )}-\frac {{St}}{\lambda _1^2}\right \vert =1. \end{eqnarray}

In figure 8(c), we plot the argument of $\chi$ as a function of $\kappa$ and ${St}$ , indicating as well the approximate bifurcation curves given by (5.14). This diagram indicates where in the $\kappa$ – ${St}$ phase space the cupula deflection follows the angular acceleration, angular velocity or the angular displacement. For small values of the Stokes number ${St}$ , we recover the previous picture, where $\kappa \ll 1$ indicates the deformation follows the angular velocity ( $\textrm{arg}(\chi )=\pi /2$ ) and $\kappa \gg 1$ indicates the response is guided by the angular acceleration ( $\textrm{arg}(\chi )=\pi$ ). However, we also find a dependence on ${St}$ for non-small Stokes numbers. For given $\kappa$ less than approximately 100, a transition to angular displacement tracking occurs for ${St}$ greater than approximately 1, meaning that the cupula system may be tuned to follow angular displacement even for small cupula stiffness if the Stokes number is high enough (characterised by $\textrm{arg}(\chi )=0$ in figure 8 c). This transition point increases for larger $\kappa$ , as indicated by the blue wedge region in figure 8(c), meaning that an orders of magnitude higher $\kappa$ is possible with cupular deflection still tracking angular velocity, if the Stokes number is accordingly increased in a very particular way. In the new regime for ${St}\gg 1$ (and moderate $\kappa$ ), we find the pressure jump is approximately

(5.15) \begin{align} \Delta p_0(t)=-\frac {\kappa }{2\alpha _0(\beta ){St}}\int _0^t\varOmega (\tau )\,\textrm {d} \tau , \end{align}

and indeed proportional to the angular displacement.

We may interpret the effect of high fluid inertia by considering the response of the cupula as the forcing frequency is increased. For small forcing frequencies, the deformation will follow the angular acceleration, and as the forcing frequency is increased, the cupula will start deforming in phase with the angular velocity, as expected from § 3.2; this is well known in the vestibular literature (Benson Reference Benson1990). However, our results from this section suggest that when the forcing frequency is further increased (in humans, to approximately 100 Hz), the cupular deformation will be in phase with angular displacement. This sort of high frequency motion could be expected for example during impacts or equipment malfunction.

5.3. Non-uniform channel widths

As noted in § 1, an advantage of our theoretical formulation is that it is also compatible with a non-uniform and arbitrary channel width, described by the function $a(s)$ , so long as the small aspect ratio between channel width and length is maintained. This case is more delicate than the one we saw in § 5.2, as the pressure gradient is no longer constant but depends on $a(s)$ , and must be integrated along the channel to obtain the pressure jump across the cupula. In Appendix D, a method based on the Laplace transform and the convolution theorem is developed, through which we obtain an approach to solving the problem for both non-uniform channel widths and ${St}\gt 0$ . The main result is that we obtain the same form of (5.7) for the deflection $\eta _0(r,t)$ , with the kernel $\mathcal{K}(x;{St})$ given by the (temporal) inverse Laplace transform of

(5.16) \begin{equation} \tilde {\mathcal{K}}(\sigma ;{St})=2\pi \left (\int _0^{2\pi }\frac {\textrm {d} s}{a(s)^4\sum _{n=1}^\infty \lambda _n^{-2}\left [a(s)^2{St}\sigma +\lambda _n^2\right ]^{-1}}\right )^{-1}. \end{equation}

For a given canal profile $a(s)$ , the Kernel may be numerically obtained by fixing a discretisation of $\sigma$ into a finite number of points. For each point, the integral can be populated for $a(s)$ and computed using standard quadrature methods. This will yield $\tilde {\mathcal{K}}(\sigma ;{St})$ for a finite number of $\sigma$ . The equation for the pressure jump (5.7) can then be either solved in Laplace space, inverting the transformed solution using an efficient algorithm (Kuhlman Reference Kuhlman2013), or in real space using a trapezoidal method.

We have developed a general framework that allows us to solve for the cupular displacement and pressure jump for complicated canal geometries allowing also for the possibility of fluid inertia. As an example of the scenarios in which this approach might be useful, we now reconsider some of the numerical results presented in figure 4.

5.4. Fluid inertia explains discrepancies between numerics and model

Figure 4 generally shows excellent agreement between the ${St}=0$ model presented in § 3 and our COMSOL numerics, especially at the scale of the largest velocities, which were used for comparison in figure 4. However, if we zoom-in to situations where the velocity is small, such as when $t=0.5$ , then we might expect to observe differences caused by small errors in the phase of the motion. Figure 9 shows just such an effect: the agreement between the predictions of the ${St}=0$ asymptotics (dashed black curves) and COMSOL simulations (blue markers) is no longer satisfactory for small and moderate values of the stiffness $\kappa \lesssim 1$ : while the absolute error is small, the relative error is very large. This is because even for small Stokes numbers, the exact time at which the velocity vanishes is different from that predicted by the ${St}=0$ asymptotics of § 3, which essentially assume that the motion is quasi-steady.

Figure 9 also shows the predictions of the leading-order in ${St}$ asymptotic results presented in this section. As might be expected, we see much better agreement between the prediction accounting for ${St}\gt 0$ via (5.2d ) and (5.7) (shown by solid red curves in figure 9) and the results of COMSOL simulations (points) than with the earlier result, which neglected the effects of fluid inertia entirely. We emphasise that this finite inertia case requires the numerical calculation of the integrals in (5.2d ) and (5.7).

Figure 9. Velocity profiles for $t=0.5$ , including the finite fluid inertia correction. Numerical results (blue markers) are compared with the theories for $St=0$ (dashed black line) and $St\gt0$ (solid red line). The parameter values are the same as in figure 7.

We now consider if domain irregularities can give rise to interesting flow phenomena, in particular, unexpected symmetry breaking and vortical flows when $\kappa \ll 1$ , which have been reported previously (Boselli et al. Reference Boselli, Obrist and Kleiser2013).

6.1. Conditions required for the formation of the utricular vortex

As noted previously, symmetry breaking in the utricle occurs when the first-order correction $\epsilon w_1$ is comparable with the leading order velocity $w_0$ . Within a cross-section, the maximum value of $w_0$ is attained at $r=0$ , and is given by

(6.2) \begin{equation} \vert w_0^*(s,t)\vert = \frac {\pi }{2I_4a(s)^2}\left \vert \dot {\varOmega }(t)+\frac {\Delta p_0}{2\pi }\right \vert \sim \frac {\pi }{2I_4a(s)^2}\vert \dot {\varOmega }(t)\vert , \end{equation}

where the last approximation follows when $\kappa \ll 1$ . The maximum value of $w_1$ is attained at $r= ({a(s)}/{\sqrt {3}})$ and $\theta =0,\pi$ , and is given by

(6.3) \begin{equation} \vert w_1^*\vert = \frac {1}{24\sqrt {3}I_4 a(s)}\big \vert -3\Delta p_0+\big(\!-6\pi +4 I_4a(s)^4 \big)\dot {\varOmega }(t)\big \vert . \end{equation}

The symmetric component of $w_1$ scales as $\epsilon \Delta p_1^{\textit{outer}} / a(s)^2$ , and since it is much smaller than both the asymmetric component of $w_1$ and the leading-order flow $w_0$ , it is neglected in the present analysis (we set $\Delta p_1^{\textit{outer}}=0$ ). Focusing on the utricle, where $a(s)$ is largest, $\vert w_1^*\vert$ is dominated by

(6.4) \begin{equation} \vert w_1^*\vert \sim \frac {a_m^3}{6\sqrt {3}}\left \vert \dot {\varOmega }(t)\right \vert . \end{equation}

Comparing $w_0$ and the next term in the expansion, $\epsilon w_1$ , we conclude that noticeably asymmetrical flow in the utricle first emerges when

(6.5) \begin{equation} \frac {\pi }{2I_4 a_m^2}\sim \epsilon \frac {a_m^3}{6\sqrt {3}},\Longleftrightarrow \frac {3\sqrt {3}\pi }{\epsilon }\sim a_m^5 \int _0^{2\pi }\frac {\textrm {d} s}{a(s)^4}. \end{equation}

Figure 11. (a) Analytical reconstruction of flow profiles in the wide region of the channel (representing the utricle), $w(r,\theta ,s=\pi ,t)$ for different values of the maximum enlargement $a_m$ . The inset shows the flow profiles in the thin region of the flow $w(r,\theta ,s=0,t)$ ; these remain symmetric, confirming the symmetry breaking mechanism is not the same as the global symmetry-breaking mechanism discussed in § 3.2. (b) Correlation (as defined in (6.6a )) between the axial velocity in the utricle $w(r)=w_0(r)+\epsilon w_1(r,\theta )$ , and the symmetric (solid) and asymmetric flow profile (dashed). Curves show the results of the analytical computation and triangles and stars show the correlations computed from the COMSOL solution. We find that the transition occurs when $\xi = a_m^5 I_4\epsilon /(3\sqrt {3}\pi ) \sim 1$ , as predicted by our analysis.

This motivates the introduction of a transition parameter $\xi = a_m^5 I_4\epsilon /(3\sqrt {3}\pi )$ . The transition can be seen qualitatively in the left panel of figure 10, where the flow in the utricle (bottom half of the torus) transitions from symmetric to asymmetric as the bulge is increased. In figure 11, this transition is analysed quantitatively. In the left panel, the flow profile in the utricle is plotted for varying $a_m$ . Observe that for large $a_m$ , the utricle flow is asymmetric, while the flow in the slender part of the canal remains symmetric. To quantify the transition to noticeably asymmetrical flow, we compute the correlation between the velocity profile $w= w_0 +\epsilon w_1$ , and the symmetric ( $w_0$ ) and asymmetric ( $w_1$ ) solutions.

The correlations may be written as

(6.6a) \begin{align} &C(w,w_0)=\frac {R(w,w_0)}{\sqrt {R(w,w)R(w_0,w_0)}},\quad C(w,w_1)=\frac {R(w,w_1)}{\sqrt {R(w,w)R(w_1,w_1)}}, \end{align}

(6.6b) \begin{align} &\qquad\qquad R(u(r,\theta ),v(r,\theta ))=\int _0^{2\pi }\int _0^{a_m}u(r,\theta )v(r,\theta )r\,\textrm {d} r\,\textrm {d} \theta , \end{align}

and we find

(6.7a) \begin{align} &\qquad\qquad C(w,w_0)=\sqrt {\frac {R_0}{R_0+\epsilon ^2 R_1}},\quad C(w,w_1)=\sqrt {\frac {\epsilon ^2 R_1}{R_0+\epsilon ^2 R_1}}, \end{align}

(6.7b) \begin{align} &\qquad R_0 = R(w_0,w_0)=2\pi \int _0^{a_m}w_0(r)^2r\,\textrm {d} r=\frac {\pi ^3}{12I_4^2a_m^2}\left (\dot {\varOmega }(t)+\frac {\Delta p_0}{2\pi }\right )^2, \end{align}

(6.7c) \begin{align} &R_1=R(w_1,w_1)=\int _0^{2\pi }\cos ^2\theta\, \textrm {d}\theta \int _0^{a_m}\left [\frac {r(r^2-a_m^2)(-3\Delta p_0+(4a_m^4I_4-6\pi )}{16I_4 a_m^4}\right ]^2r\,\textrm {d} r \nonumber \\ &\quad \,=\frac {\pi }{24\boldsymbol{\cdot }16^2I_4^2}\big (-3\Delta p_0+(4a_m^4I_4-6\pi )\dot {\varOmega }(t)\big )^2. \end{align}

Here, we have assumed that $\Delta p_1^{\textit{outer}}=0$ and used the fact that $R(w_0,w_1)=0$ , which is trivial since $w_0w_1\sim \cos \theta$ . The correlation $C(w,w_0)$ will be close to 1 when the flow is largely symmetrical, and closer to zero when the flow is noticeably asymmetrical. The opposite is true of the correlation $C(w,w_1)$ . In figure 11(a), we plot the symmetrical correlation $C(w,w_0)$ (solid line) and the asymmetrical correlation $C(w,w_1)$ (dashed line) as predicted from our analytical model.

Based on the definition of $\xi$ , we expect a transition when $\xi \sim 1$ . Figure 11 confirms this expectation, showing that a transition indeed occurs at this point when the geometry $a(s)$ satisfies (6.1). We note that for times when $\dot {\varOmega }(t) = 0$ , our analysis does not apply and the flow remains symmetric, even for large utricles. This follows from (6.4), as it is clear that $w_1 = 0$ when $\dot {\varOmega }(t) = 0$ .

Our analysis of the onset of vortical flow is further supported by numerical simulations, represented by triangles and stars in figure 11(b). These simulations were conducted using COMSOL, modelling the utricle as an ellipsoidal expansion of the toroidal geometry. We find good agreement between the analytical predictions and COMSOL simulations regarding the emergence of vortical flow. (We attribute the small discrepancy between the COMSOL and analytical model to the fact that the geometry is not identical in both cases, and furthermore, because the analytical result does not have a solid obstacle disrupting the flow.) Furthermore, the consistency between the numerical simulations and analytical solutions suggests that, although the first-order correction may be larger than the leading-order term, the subsequent terms in the series remain well behaved, ensuring that the asymptotic ordering is preserved.