3.1. The duct outlet as an acoustic source

Here, we wish to model the propagation of sound from a duct outlet, in order to inform our choice of radiation condition for the duct itself. Rather than a duct radiating sound into free space, we consider an inner duct positioned within an outer duct several times larger, i.e. ‘approximately free’ space. A schematic of this model in depicted in figure 1. The inner-duct exit, shown in green in figure 1, is modelled as an acoustic source within the outer duct at arclength $s = 0$ , with its width being a fraction $\eta \lt 1$ of the outer duct’s width. The right (positive $s$ -facing) side of the source then emanates sound into the outer duct, while the left side (shown in red in figure 1) will have an acoustic absorption condition, to be specified later: since sound entering this surface departs forever from the domain of computation, we term this the ‘acoustic wormhole’. This is obviously a somewhat idealised situation, as the body of the inner duct is not present in the outer domain, but the expectation is as follows: of the sound leaving the bell of a trumpet, only a relatively small proportion will emanate backwards; a yet smaller proportion of this will reflect off the body of the trumpet; and this small fraction of reflected sound’s contribution is unlikely to have a significant influence on the inner duct’s acoustics. While there are many possible boundary conditions that could be applied to the outer duct, here, we will use a rigid wall boundary condition, as is used for the inner duct. This assumption will be justified and discussed further in § 5.1.3. As an idea of the physical parameters typically involved, a brass instrument has a bell that is typically around 10 cm or so in diameter. That means, with an $\eta = 1/40$ ratio (used for end-correction calculations later on) the outer duct is 4 m in diameter. Clearly, unless the player is outside, reflections from the room the player is in would be comparable to the reflections from the ‘fake’ outer duct employed here, and so are unlikely to affect the inner-duct resonances. The curved arrows in figure 1 demonstrate that sound (i) leaves the inner-duct domain via the green surface, (ii) passes into the right-hand side of the outer duct, (iii) can also enter the left-hand side of the outer duct and (iv) may be absorbed by the ‘acoustic wormhole’ surface shown in red.

We define $p^R(\boldsymbol{x},t)$ and $u^R(\boldsymbol{x},t)$ to be the pressure and longitudinal velocity just to the right of the source ( $s=0^+$ ), while $p^L(\boldsymbol{x},t)$ and $u^L(\boldsymbol{x},t)$ are the equivalent quantities just to the left of the source ( $s=0^-)$ . These quantities are all expanded in the same way as described in § 2, although being within the outer duct they are all expanded in terms of the outer-duct modes: taking $p^R(\boldsymbol{x},t)$ as an example, we have

(3.1) \begin{equation} p^R(\boldsymbol{x},t) = \sum _{a = -\infty }^\infty \sum _{\alpha = 0}^{\infty }{\mathrm{e}}^{-{i} a\omega t}\psi _\alpha ^2(\boldsymbol{x})P_{\alpha ,2}^{\,a,R}, \end{equation}

where $\psi _\alpha ^2(\boldsymbol{x})$ are the outer-duct modes given by (2.3) in two dimensions and (2.4) in three dimensions (using the outer-duct width and $s = 0$ in these equations). Note that the definition of the coefficient $P_{\alpha ,2}^{\,a,R}$ is

(3.2) \begin{equation} P_{\alpha ,2}^{\,a,R} = \frac {\omega }{2\pi }\int _{t = 0}^{2\pi }\iint _{S_2}p^R(\boldsymbol{x},t){\mathrm{e}}^{{i} a\omega t}\psi _\alpha ^2(\boldsymbol{x})\mathrm{d}S\mathrm{d}t, \end{equation}

where $S_2$ is the outer-duct cross-section at 0.

We then treat $P_{\alpha ,2}^{\,a,R}$ as a set of vectors with $\alpha$ -coefficients, $\{\boldsymbol{p}_2^{\,a,R}:\, a\in \mathbb{Z}\}$ , which we refer to with the shorthand notation $\boldsymbol{p}_2^R$ (and similarly with $\boldsymbol{u}_2^R$ , $\boldsymbol{p}_2^L$ and $\boldsymbol{u}_2^L$ , an equivalent shorthand to the one seen in (2.6)). We can then introduce tensorial linear and nonlinear admittances to the right of the source $(s = 0^+)$ , ${\unicode{x1D654}}_{\alpha \beta }^{\,a,R}$ and $\mathcal{Y}_{\alpha \beta \gamma }^{\,ab,R}$ , defined as in (2.6), and similarly for $s = 0^-$ , giving ${\unicode{x1D654}}_{\alpha \beta }^{\,a,L}$ and $\mathcal{Y}_{\alpha \beta \gamma }^{\,ab,L}$ . Equipped with these radiative quantities, we then stipulate that sound to the right of the source must only travel rightward in the outer duct, and sound to the left must only travel leftward; equivalently, we are taking the outer-duct admittances ( ${\unicode{x1D654}}^{\,R}_2$ , $\mathcal{Y}^R_2$ ) to be the positive characteristic ( $\overline {{\unicode{x1D654}}}_2$ , $\overline {\mathcal{Y}}_2$ ) for $s \gt 0$ , and ${\unicode{x1D654}}^{\,L}_2$ , $\mathcal{Y}^L_2$ to be the negative characteristic ( $-\overline {{\unicode{x1D654}}}_2$ , $-\overline {\mathcal{Y}}_2$ ) for $s \lt 0$ . Note that the impedances in these regions will then be ( $\overline {{\unicode{x1D655}}}_2,\overline {\mathcal{Z}}_2$ ) and ( $-\overline {{\unicode{x1D655}}}_2,\overline {\mathcal{Z}}_2$ ) respectively (the non-matching signs here are due to reflection in $s=0$ giving a sign change of longitudinal velocity $u'$ but not of pressure $p'$ , and the linear and weakly nonlinear parts of the impedance being either linear or quadratic in this sign change), where $\overline {{\unicode{x1D655}}}_2 = {\overline {{\unicode{x1D654}}}_2}^{-1}$ and $\overline {\mathcal{Z}}_2 = -\overline {{\unicode{x1D655}}}_2\overline {\mathcal{Y}}_2\langle \overline {{\unicode{x1D655}}}_2,\overline {{\unicode{x1D655}}}_2\rangle$ . This notation is taken from Jensen & Brambley (Reference Jensen and Brambley2026), where explicit formulae and a derivation of the characteristic admittance is given.

3.2. Restriction and zero extension

Before we derive the outlet condition, it will first prove helpful to consider how modal expansions using the outer-duct modes $\psi _\alpha ^2(\boldsymbol{x})$ are related to modal expansions using the inner-duct modes; inner-duct modes are denoted here as $\psi _\alpha ^1(\boldsymbol{x})$ , and are similarly given by (2.3) using the inner-duct width $X_1$ in two dimensions and by (2.4) using the inner-duct radius $R_1$ in three dimensions. Expanding using the inner-duct modes would result in a different set of coefficients $\boldsymbol{p}_1^R$ , $\boldsymbol{u}_1^R$ , $\boldsymbol{p}_1^L$ and $\boldsymbol{u}_1^L$ , and is valid only for $|x| \lt X_1/2$ in two dimensions and $r \lt R_1$ in three dimensions. We now relate the outer-duct coefficients to the inner-duct ones. For example,

(3.3) \begin{equation} \sum _{\beta = 0}^\infty \psi _\beta ^1(\boldsymbol{x})P_{\beta ,1}^{\,a,R} = \sum _{\beta = 0}^\infty \psi _\beta ^2(\boldsymbol{x})P_{\beta ,2}^{\,a,R} \qquad \text{provided}\quad \boldsymbol{x}\in S_1, \end{equation}

where $S_1$ is the inner-duct cross-section at $s = 0$ . Multiplying both sides by $\psi _\alpha ^1(\boldsymbol{x})$ and integrating both sides over an inner-duct cross-section gives (by orthogonality of the modes)

(3.4) \begin{equation} P_{\alpha ,1}^{\,a,R} = \sum _{\beta = 0}^\infty \int _{S_1}\psi _\alpha ^1(\boldsymbol{x})\psi _\beta ^2(\boldsymbol{x})\mathrm{d}S\,P_{\beta ,2}^{\,a,R} =: \sum _\beta {\unicode{x1D641}}_{\alpha \beta } P_{\beta ,2}^{\,a,R}, \end{equation}

where we term $\unicode{x1D641}$ the restriction operator. We utilise the vector notation here and write $\boldsymbol{p}_1^R = {\unicode{x1D641}}\boldsymbol{p}_2^R$ .

Since restriction is a one-way operation (in that information for $\boldsymbol{x}\in S_2\setminus S_1$ is lost), we cannot define an inverse operation, but we can still map from the inner duct to the outer duct if we think carefully about what this lost information can be replaced with. While we expect a jump in the pressure and velocity as we move from one side of the source to the other, we do not expect this discontinuity to hold around the source. We therefore define acoustic jump variables $u^*(\boldsymbol{x},t)$ and $p^*(\boldsymbol{x},t)$ as

(3.5) \begin{align} u^*(\boldsymbol{x},t) = u^R(\boldsymbol{x},t)|_{s = 0^+} - u^L(\boldsymbol{x},t)|_{s = 0^-}, && p^*(\boldsymbol{x},t) = p^R(\boldsymbol{x},t)|_{s = 0^+} - p^L(\boldsymbol{x},t)|_{s = 0^-}, \end{align}

in the expectation that these will be non-zero only on the region $S_1$ . Finding the outer-duct coefficients of one of these jump variables, we have

(3.6) \begin{equation} U_{\alpha ,2}^{a,*} = \frac {\omega }{2\pi }\int _{t = 0}^{2\pi }\iint _{S_2}u^*(\boldsymbol{x},t){\mathrm{e}}^{{i} a\omega t}\psi _\alpha ^2(\boldsymbol{x})\mathrm{d}S\mathrm{d}t, \end{equation}

but the integral over $S_2$ here may equivalently be performed over $S_1$ , since $u^*$ is zero outside $S_1$ . Making this change, and expanding $u^*(\boldsymbol{x},t)$ in terms of the inner-duct modes, we get

(3.7) \begin{equation} U_{\alpha ,2}^{a,*} = \iint _{S_1}\psi _\alpha ^2(\boldsymbol{x})\psi _\beta ^1(\boldsymbol{x})\mathrm{d}S\,U_{\beta ,1}^{a,*}, \end{equation}

or in other words

(3.8) \begin{align} \boldsymbol{u}_2^* = {\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^*, &&\text{and by the same property,} &&\boldsymbol{p}_2^* = {\unicode{x1D641}}^{{T}}\!\boldsymbol{p}_1^*. \end{align}

From this, we deduce that taking the transpose of the restriction operator results in a zero-extension operator. Note that, while ${\unicode{x1D641}}{\unicode{x1D641}}^{{T}}\! = {\unicode{x1D644}}$ (because zero extending and then restricting returns the original), in general ${\unicode{x1D641}}^{{T}}\!{\unicode{x1D641}}\! \neq {\unicode{x1D644}}$ (since restricting and then zero extending essentially sets the original to zero for $\boldsymbol{x}\in S_2\setminus S_1$ ).

For the mode definitions (2.3) and (2.4), the restriction operator is given in two dimensions by

(3.9) \begin{equation} {\unicode{x1D641}}_{\alpha \beta } = \begin{cases} \dfrac {C_\alpha C_\beta }{\pi }\dfrac {\eta ^{3/2}\beta }{\alpha ^2 - \eta ^2\beta ^2}\left [\sin \left (\dfrac {1 - \eta }{2}\beta \pi \right ) - (-1)^\alpha \sin \left (\dfrac {1 + \eta }{2}\beta \pi \right )\right ] &\text{if }\alpha \neq \eta \beta , \\[15pt] \dfrac {C_\alpha C_\beta }{2}\sqrt {\eta }\cos \left (\dfrac {\alpha \pi }{2}\left [1 - \dfrac {1}{\eta }\right ]\right ) &\text{if }\alpha = \eta \beta \neq 0, \\[15pt] C_\alpha C_\beta \sqrt {\eta } &\text{if }\alpha = \eta \beta = 0, \end{cases} \end{equation}

and in three dimensions by

(3.10) \begin{equation} {\unicode{x1D641}}_{\alpha \beta } = \begin{cases} (1+\delta _{m_\alpha 0})\eta C_\alpha C_\beta &\dfrac {\lambda _\alpha \mathrm{J}_{m_\alpha - 1}(\lambda _\alpha )\mathrm{J}_{m_\alpha }(\eta \lambda _\beta ) - \eta \lambda _\beta \mathrm{J}_{m_\alpha }(\lambda _\alpha )\mathrm{J}_{m_\alpha - 1}(\eta \lambda _\beta )}{\eta ^2\lambda _\beta ^2 - \lambda _\alpha ^2} \\[15pt] &\qquad \qquad \qquad \qquad \qquad \text{if }m_\alpha = m_\beta ,\,\xi _\alpha = \xi _\beta ,\,\eta \lambda _\beta \neq \lambda _\alpha ,\\[9pt] \eta & \qquad \qquad \qquad \qquad \qquad \text{if }m_\alpha = m_\beta ,\,\xi _\alpha = \xi _\beta ,\,\eta \lambda _\beta = \lambda _\alpha , \\[9pt] 0 & \qquad \qquad \qquad \qquad \qquad \text{otherwise}. \end{cases} \end{equation}

3.3. Derivation of the outlet admittance

We now go on to derive the open-end admittance boundary condition, using the restriction and zero-extension operators $\unicode{x1D641}$ and ${\unicode{x1D641}}^{{T}}$ from the previous section, applied at $s = 0$ . Specifically, we aim to derive a condition relating pressure and axial velocity at the exit of the smaller duct as it enters the larger duct (shown in green in figures 1 a and 1 b); i.e. a relation between $\boldsymbol{p}_1^R$ and $\boldsymbol{u}_1^R$ at $s=0$ . To do so, we will assume: (i) a known relation between $\boldsymbol{p}_1^L$ and $\boldsymbol{u}_1^L$ (i.e. the condition on the red ‘acoustic wormhole’ surface in figure 1 b); (ii) that waves are only outgoing in the larger duct for both $s\gt 0$ and $s\lt 0$ ; and (iii) that the acoustic variables are continuous at the duct exit (across $s=0$ ) in the larger duct outside the smaller duct (for $\boldsymbol{x}\in S_2\setminus S_1$ ). This last condition means that the jump in velocity and pressure when crossing $s=0$ , given by $u^*$ and $p^*$ in (3.5), is only non-zero for $\boldsymbol{x}\in S_1$ . The entirety of this section is evaluated within the outer duct at $s=0$ . We start by using the acoustic jump variables in the outer duct to eliminate $\boldsymbol{u}_2^L = \boldsymbol{u}_2^R - \boldsymbol{u}_2^*$ and $\boldsymbol{p}_2^L = \boldsymbol{p}_2^R - \boldsymbol{p}_2^*$ ; in particular, the outer-left-duct admittance equation becomes

(3.11a) \begin{align} \boldsymbol{u}_2^R - \boldsymbol{u}_2^* &= -\overline {{\unicode{x1D654}}}_2\big (\boldsymbol{p}_2^R - \boldsymbol{p}_2^*\big ) - \overline {\mathcal{Y}}_2\big \langle \boldsymbol{p}_2^R - \boldsymbol{p}_2^*,\boldsymbol{p}_2^R - \boldsymbol{p}_2^*\big \rangle , \end{align}

or equivalently,

(3.11b) \begin{align} \boldsymbol{p}_2^R - \boldsymbol{p}_2^* &= -\overline {{\unicode{x1D655}}}_2\big (\boldsymbol{u}_2^R - \boldsymbol{u}_2^*\big ) + \overline {\mathcal{Z}}_2\big \langle \boldsymbol{u}_2^R - \boldsymbol{u}_2^*,\boldsymbol{u}_2^R - \boldsymbol{u}_2^*\big \rangle . \end{align}

At linear order, this may be rearranged for $\boldsymbol{u}_2^R$ using the identities $\boldsymbol{u}_2^R = \overline {{\unicode{x1D654}}}_2\boldsymbol{p}_2^R$ and $\boldsymbol{p}_2^R = \overline {{\unicode{x1D655}}}_2\boldsymbol{u}_2^R$ to give

(3.12) \begin{align} &\boldsymbol{u}_2^R = \frac {1}{2} \big (\boldsymbol{u}_2^* + \overline {{\unicode{x1D654}}}_2\boldsymbol{p}_2^*\big ) +O(M^2),&&\text{or equivalently}&&\boldsymbol{p}_2^R = \frac {1}{2} \big (\boldsymbol{p}_2^* + \overline {{\unicode{x1D655}}}_2\boldsymbol{u}_2^*\big )+O(M^2). \end{align}

These two expressions may be substituted into the nonlinear terms in (3.11) (recall we are neglecting $O(M^3)$ and smaller terms), resulting in

(3.13a) \begin{align} \boldsymbol{u}_2^R &= \frac {1}{2} \big (\boldsymbol{u}_2^* + \overline {{\unicode{x1D654}}}_2\boldsymbol{p}_2^*\big ) + \frac {1}{8}\overline {\mathcal{Y}}_2\Big (\big \langle \boldsymbol{p}_2^* + \overline {{\unicode{x1D655}}}_2\boldsymbol{u}_2^*,\boldsymbol{p}_2^* + \overline {{\unicode{x1D655}}}_2\boldsymbol{u}_2^*\big \rangle - \big \langle \boldsymbol{p}_2^* - \overline {{\unicode{x1D655}}}_2\boldsymbol{u}_2^*,\boldsymbol{p}_2^* - \overline {{\unicode{x1D655}}}_2\boldsymbol{u}_2^*\big \rangle \Big ), \\[-12pt]\nonumber \end{align}

(3.13b) \begin{align} \boldsymbol{p}_2^R &= \frac {1}{2} \big (\boldsymbol{p}_2^* + \overline {{\unicode{x1D655}}}_2\boldsymbol{u}_2^*\big ) + \frac {1}{8}\overline {\mathcal{Z}}_2\Big (\big \langle \boldsymbol{u}_2^* + \overline {{\unicode{x1D654}}}_2\boldsymbol{p}_2^*,\boldsymbol{u}_2^* + \overline {{\unicode{x1D654}}}_2\boldsymbol{p}_2^*\big \rangle + \big \langle \boldsymbol{u}_2^* - \overline {{\unicode{x1D654}}}_2\boldsymbol{p}_2^*,\boldsymbol{u}_2^* - \overline {{\unicode{x1D654}}}_2\boldsymbol{p}_2^*\big \rangle \Big ). \end{align}

By outer-duct continuity at $s=0$ , both $u^*$ and $p^*$ are zero outside $\boldsymbol{x}\in S_1$ , and so we can use (3.8) to replace all outer-duct quantities on the right-hand sides with inner-duct ones. Applying $\unicode{x1D641}$ to both sides then gives

(3.14a) \begin{align} \boldsymbol{u}_1^R = &\frac {1}{2}\big (\boldsymbol{u}_{1}^* + {\unicode{x1D641}}\overline {{\unicode{x1D654}}}_2{\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_{1}^*\big ) \\ \notag &+ \frac {1}{8}{\unicode{x1D641}}\overline {\mathcal{Y}}_2\Big (\big \langle {\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_1^* + \overline {{\unicode{x1D655}}}_2{\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^*,{\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_1^* + \overline {{\unicode{x1D655}}}_2{\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^*\big \rangle - \big \langle {\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_1^* - \overline {{\unicode{x1D655}}}_2{\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^*,{\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_1^* - \overline {{\unicode{x1D655}}}_2{\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^*\big \rangle \Big ), \\[-12pt]\nonumber \end{align}

(3.14b) \begin{align} \boldsymbol{p}_1^R = &\frac {1}{2}\big (\boldsymbol{p}_1^* + {\unicode{x1D641}}\overline {{\unicode{x1D655}}}_2{\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^*\big ) \\ \notag &+ \frac {1}{8}{\unicode{x1D641}}\overline {\mathcal{Z}}_2\Big (\big \langle {\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^* + \overline {{\unicode{x1D654}}}_2{\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_1^*,{\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^* + \overline {{\unicode{x1D654}}}_2{\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_1^*\big \rangle + \big \langle {\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_{1}^* - \overline {{\unicode{x1D654}}}_2{\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_1^*,{\unicode{x1D641}}^{{T}}\!\boldsymbol{u}_1^* - \overline {{\unicode{x1D654}}}_2{\unicode{x1D641}}^{{T}}\!\!\boldsymbol{p}_1^*\big \rangle \Big ). \end{align}

The inner-duct admittances to the right of the source, ${\unicode{x1D654}}_1^{\,R}$ and $\mathcal{Y}_1^R$ , are what we wish to determine; those to the left of the source, ${\unicode{x1D654}}_1^{\,L}$ and $\mathcal{Y}_1^L$ , will be prescribed. For the remainder of this section, we will suppress indices on radiative quantities and (unless otherwise specified) take $({\unicode{x1D654}}^{\,R},\mathcal{Y}^R)$ to be $({\unicode{x1D654}}_1^{\,R},\mathcal{Y}_1^R)$ , $({\unicode{x1D654}}^{\,L},\mathcal{Y}^L)$ to be $({\unicode{x1D654}}_1^{\,L},\mathcal{Y}_1^L)$ and $(\overline {{\unicode{x1D654}}},\overline {\mathcal{Y}})$ to be $(\overline {{\unicode{x1D654}}}_2,\overline {\mathcal{Y}}_2)$ (with similar shorthand for all impedances). We then have the acoustic wormhole’s impedance equation in terms of the jump and outlet quantities

(3.15) \begin{equation} \boldsymbol{p}_1^* - \boldsymbol{p}_1^R = {\unicode{x1D655}}^L\big (\boldsymbol{u}_1^* - \boldsymbol{u}_1^R\big ) - \mathcal{Z}^L\big \langle \boldsymbol{u}_1^* - \boldsymbol{u}_1^R,\boldsymbol{u}_1^* - \boldsymbol{u}_1^R\big \rangle . \end{equation}

Substituting $\boldsymbol{p}_1^*$ from this expression into (3.14a ), at linear order we may rearrange to find

(3.16) \begin{equation} \boldsymbol{u}_1^* = ({\unicode{x1D644}} + {\unicode{x1D64C}})\boldsymbol{u}_1^R - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R + O\big (M^2\big ), \end{equation}

where ${\unicode{x1D64C}} := ({\unicode{x1D644}} + {\unicode{x1D641}}\overline {{\unicode{x1D654}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D655}}^L )^{-1}$ is a shorthand quantity. This expression can be substituted back into (3.15) to gain an expression for $\boldsymbol{p}_1^*$ , again at linear order

(3.17) \begin{equation} \boldsymbol{p}_1^* = {\unicode{x1D655}}^L{\unicode{x1D64C}}\big (\boldsymbol{u}_1^R + {\unicode{x1D654}}^{\,L}\!\boldsymbol{p}_1^R\big ) + O\big (M^2\big ). \end{equation}

We may now perform the same two steps, but including the quadratic quantities, making use of these linear identities within the nonlinear terms (since the error introduced by doing so is $O (M^3 )$ , i.e. comparable to the magnitude of other neglected terms). This results in

(3.18a) \begin{align} \boldsymbol{u}_1^* = ({\unicode{x1D644}} + {\unicode{x1D64C}})\boldsymbol{u}_1^R - &({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R + ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\mathcal{Z}^L\big \langle {\unicode{x1D64C}}\boldsymbol{u}_1^R - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R,{\cdots} \big \rangle \\ \notag + \frac {1}{4}{\unicode{x1D64C}}{\unicode{x1D641}}\overline {{\unicode{x1D654}}}\overline {\mathcal{Z}}\big (- &\,\big \langle \overline {{\unicode{x1D654}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D655}}^L{\unicode{x1D64C}}\big (\boldsymbol{u}_1^R + {\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R\big ) - {\unicode{x1D641}}^{{T}}\big [({\unicode{x1D644}} + {\unicode{x1D64C}})\boldsymbol{u}_1^R - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R\big ],{\cdots} \big \rangle \\ \notag +&\,\big \langle \overline {{\unicode{x1D654}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D655}}^L{\unicode{x1D64C}}\big (\boldsymbol{u}_1^R + {\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R\big ) + {\unicode{x1D641}}^{{T}}\big [({\unicode{x1D644}} + {\unicode{x1D64C}})\boldsymbol{u}_1^R - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R\big ],{\cdots} \big \rangle \big ), \end{align}

(3.18b) \begin{align} \boldsymbol{p}_1^* &= {\unicode{x1D655}}^L{\unicode{x1D64C}}\big (\boldsymbol{u}_1^R + {\unicode{x1D654}}^{\,L} \boldsymbol{p}_1^R\big ) - {\unicode{x1D655}}^L{\unicode{x1D64C}}{\unicode{x1D654}}^{\,L}\mathcal{Z}^L\big \langle {\unicode{x1D64C}}\boldsymbol{u}_1^R - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R,{\cdots} \big \rangle \\ \notag &\,+ \frac {1}{4}{\unicode{x1D655}}^L{\unicode{x1D64C}}{\unicode{x1D641}}\overline {{\unicode{x1D654}}}\overline {\mathcal{Z}}\big (- \big \langle \overline {{\unicode{x1D654}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D655}}^L{\unicode{x1D64C}}\big (\boldsymbol{u}_1^R + {\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R\big ) - {\unicode{x1D641}}^{{T}}\big [({\unicode{x1D644}} + {\unicode{x1D64C}})\boldsymbol{u}_1^R - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R\big ],{\cdots} \big \rangle \\ \notag &\qquad \qquad \quad + \big \langle \overline {{\unicode{x1D654}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D655}}^L{\unicode{x1D64C}}\big (\boldsymbol{u}_1^R + {\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R\big ) + {\unicode{x1D641}}^{{T}}\big [({\unicode{x1D644}} + {\unicode{x1D64C}})\boldsymbol{u}_1^R - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R\big ],{\cdots} \big \rangle \big ), \end{align}

where ‘ ${\cdots}$ ’ indicates the left argument of the $\langle \bullet ,\bullet \rangle$ operation being repeated for the right. Since we have not yet used (3.14b ), we may now do so, and so we substitute (3.18b ) and (3.18a ) into it. At linear order, this results in

(3.19) \begin{align} \boldsymbol{p}_1^R = \frac {1}{2}\big [{\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!({\unicode{x1D644}} + {\unicode{x1D64C}}) + {\unicode{x1D655}}^L{\unicode{x1D64C}}\big ]\boldsymbol{u}_1^R + \frac {1}{2}\big [-{\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!({\unicode{x1D644}} - {\unicode{x1D64C}}) + {\unicode{x1D655}}^L{\unicode{x1D64C}}\big ]{\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^R + O\big (M^2\big ). \end{align}

This may be rearranged, using another shorthand quantity $\widetilde {{\unicode{x1D64C}}} := ({\unicode{x1D644}} + {\unicode{x1D654}}^{\,L}{\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}} )^{-1}$ , into the form $\boldsymbol{p}_1^R$ = ${\unicode{x1D655}}^R\boldsymbol{u}_1^R + O (M^2 )$ , where

(3.20) \begin{equation} {\unicode{x1D655}}^R = -{\unicode{x1D655}}^L\big [{\unicode{x1D644}} + \widetilde {{\unicode{x1D64C}}}^{-1}({\unicode{x1D644}} - {\unicode{x1D64C}})\big ]^{-1}\big [{\unicode{x1D644}} - \widetilde {{\unicode{x1D64C}}}^{-1}({\unicode{x1D644}} + {\unicode{x1D64C}})\big ]. \end{equation}

As before, we now carry out the same steps, but including the quadratic quantities, and making use of linear identities within the nonlinear terms. The full nonlinear substitution of (3.18b ) and (3.18a ) into (3.14b ) is then $\boldsymbol{p}_1^R = {\unicode{x1D655}}^R\boldsymbol{u}_1^R + \mathcal{Z}^R\langle \boldsymbol{u}_1^R,\boldsymbol{u}_1^R\rangle$ , where

(3.21) \begin{align} \mathcal{Z}^R = {\unicode{x1D655}}^L\big [{\unicode{x1D644}} + \widetilde {{\unicode{x1D64C}}}^{-1}({\unicode{x1D644}} - {\unicode{x1D64C}})\big ]^{-1}&\bigg \{\!-\big [{\unicode{x1D644}} - \widetilde {{\unicode{x1D64C}}}^{-1}({\unicode{x1D644}} - {\unicode{x1D64C}})\big ]{\unicode{x1D654}}^{\,L}\mathcal{Z}^L\Big \langle {\unicode{x1D64C}} - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}{\unicode{x1D655}}^R,{\cdots} \!\Big \rangle \\ \notag + \frac {1}{4}\big (\!{\unicode{x1D654}}^{\,L}{\unicode{x1D641}} - \widetilde {{\unicode{x1D64C}}}^{-1}{\unicode{x1D64C}}{\unicode{x1D641}}\overline {{\unicode{x1D654}}}\big )\, \overline {\mathcal{Z}}\, &\Big \langle \overline {{\unicode{x1D654}}}{\unicode{x1D641}}^{{T}}{\unicode{x1D655}}^L{\unicode{x1D64C}}\big ({\unicode{x1D644}} + \!{\unicode{x1D654}}^{\,L}{\unicode{x1D655}}^R\big ) - {\unicode{x1D641}}^{{T}}\big [({\unicode{x1D644}} +\! {\unicode{x1D64C}}) - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}{\unicode{x1D655}}^R\big ],{\cdots} \!\Big \rangle \\ \notag + \frac {1}{4}\big (\!{\unicode{x1D654}}^{\,L}{\unicode{x1D641}} + \widetilde {{\unicode{x1D64C}}}^{-1}{\unicode{x1D64C}}{\unicode{x1D641}}\overline {{\unicode{x1D654}}}\big )\, \overline {\mathcal{Z}} \, &\Big \langle \overline {{\unicode{x1D654}}}{\unicode{x1D641}}^{{T}}{\unicode{x1D655}}^L{\unicode{x1D64C}}\big ({\unicode{x1D644}} + \!{\unicode{x1D654}}^{\,L}{\unicode{x1D655}}^R\big ) + {\unicode{x1D641}}^{{T}}\big [({\unicode{x1D644}} +\! {\unicode{x1D64C}}) - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}{\unicode{x1D655}}^R\big ],{\cdots} \!\Big \rangle \bigg \}, \end{align}

and we recall the notation that $\mathcal{A} = \mathcal{B}\langle {\unicode{x1D63E}},{\unicode{x1D63F}}\rangle$ defines a nonlinear operator such that $\mathcal{A}\langle \boldsymbol{v}, \boldsymbol{w}\rangle = \mathcal{B}\langle {\unicode{x1D63E}}\boldsymbol{v}, {\unicode{x1D63F}}\boldsymbol{w}\rangle$ . Using the nonlinear inversion rule from § 2, we may invert the impedances ${\unicode{x1D655}}^R$ and $\mathcal{Z}^R$ to get the admittances

(3.22a) \begin{align} & \hspace{-5.5cm} {\unicode{x1D654}}^{\,R} = -\big [{\unicode{x1D644}} - \widetilde {{\unicode{x1D64C}}}^{-1}({\unicode{x1D644}} + {\unicode{x1D64C}})\big ]^{-1}\big [{\unicode{x1D644}} + \widetilde {{\unicode{x1D64C}}}^{-1}({\unicode{x1D644}} - {\unicode{x1D64C}})\big ]{\unicode{x1D654}}^{\,L}, \end{align}

(3.22b) \begin{align} \mathcal{Y}^R =& \big [{\unicode{x1D644}} - \widetilde {{\unicode{x1D64C}}}^{-1}({\unicode{x1D644}} + {\unicode{x1D64C}})\big ]^{-1}\bigg \{\big [{\unicode{x1D644}} - \widetilde {{\unicode{x1D64C}}}^{-1}({\unicode{x1D644}} - {\unicode{x1D64C}})\big ]\mathcal{Y}^L\Big \langle {\unicode{x1D644}} - {\unicode{x1D655}}^L{\unicode{x1D64C}}\left ({\unicode{x1D654}}^{\,L} + {\unicode{x1D654}}^{\,R}\right ),{\cdots} \Big \rangle \\ \notag & - \frac {1}{4}\big (\!{\unicode{x1D654}}^{\,L}{\unicode{x1D641}}\overline {{\unicode{x1D655}}} - \widetilde {{\unicode{x1D64C}}}^{-1}{\unicode{x1D64C}}{\unicode{x1D641}}\big )\overline {\mathcal{Y}}\Big \langle {\unicode{x1D641}}^{{T}}{\unicode{x1D655}}^L{\unicode{x1D64C}}\big ({\unicode{x1D654}}^{\,L} + \!{\unicode{x1D654}}^{\,R}\big ) - \overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\big [({\unicode{x1D644}} + {\unicode{x1D64C}}){\unicode{x1D654}}^{\,R} - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\big ],{\cdots} \Big \rangle \\ \notag & - \frac {1}{4}\big (\!{\unicode{x1D654}}^{\,L}{\unicode{x1D641}}\overline {{\unicode{x1D655}}} + \widetilde {{\unicode{x1D64C}}}^{-1}{\unicode{x1D64C}}{\unicode{x1D641}}\big )\overline {\mathcal{Y}}\Big \langle {\unicode{x1D641}}^{{T}}{\unicode{x1D655}}^L{\unicode{x1D64C}}\big ({\unicode{x1D654}}^{\,L} + \!{\unicode{x1D654}}^{\,R}\big ) + \overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\big [({\unicode{x1D644}} + {\unicode{x1D64C}}){\unicode{x1D654}}^{\,R} - ({\unicode{x1D644}} - {\unicode{x1D64C}}){\unicode{x1D654}}^{\,L}\big ],{\cdots} \Big \rangle \bigg \}. \end{align}

It finally remains to choose the boundary condition to apply on the acoustic wormhole (left of the discontinuity) by specifying the admittance ${\unicode{x1D654}}^{\,L}$ and $\mathcal{Y}^L$ . This boundary is somewhat artificial, as we have here neglected the outer walls of the inner duct, and this boundary sits within the region that should be inaccessible to sound from outside the inner duct (see figure 1 b for a corresponding schematic). In order to avoid spurious reflections from this artificial boundary, we here choose the left-hand side admittances to be perfectly absorbing, i.e. ${\unicode{x1D654}}^{\,L} = \overline {{\unicode{x1D654}}}_1$ and $\mathcal{Y}^L = \overline {\mathcal{Y}}_1$ , and hence sound may enter the ‘acoustic wormhole’ but no sound leaves from it. This, in combination with (3.22a ) and (3.22b ), fully determines the admittance at the inner duct’s open end, which may be used as an initial condition to solve the Riccati-style system in $({\unicode{x1D654}}_1(s),\mathcal{Y}_1(s))$ within the inner duct from arclength $s_{\text{o}}$ to $0$ (as described in Jensen & Brambley Reference Jensen and Brambley2026), so that the admittance is then known everywhere throughout the system, and can then be used to calculate acoustic variables such as the pressure.

3.4. Pressure continuity conditions

Once the admittance has been calculated in the inner duct, the pressure in the inner duct $\boldsymbol{p}_1^R(s)$ may be computed (subject to an inlet condition). Details of this calculation are the subject of Jensen & Brambley (Reference Jensen and Brambley2026). For the examples that follow, we will specify the inlet condition in terms of the incoming (forward-going) acoustics only; that is, we assume that acoustics travelling backwards along the duct away from the open end may leave the duct at the inlet unimpeded. While this is a helpful inlet condition to use for studying duct resonances, as we will do later, we note in passing that it is equally possible to use this framework for any inlet condition, or indeed for a duct open at both ends with a suitable sound source located in its interior (which would involve the decomposition of the problem into two ducts either side of the source, with the two sides potentially communicating with one another across it).

In this section, we derive the relations needed for transforming the computed pressure at the outlet of the inner duct into a boundary condition for computation of the pressure field in the outer duct. We start by imposing continuity of acoustic velocity at $s=0$ outside the smaller duct, so that the velocity discontinuity is zero extended from the outlet

(3.23) \begin{equation} \boldsymbol{u}_2^R - \boldsymbol{u}_2^L = {\unicode{x1D641}}^{{T}}\!\big (\boldsymbol{u}_1^R - \boldsymbol{u}_1^L\big ). \end{equation}

In terms of the pressure, using the outer-duct admittance gives

(3.24) \begin{align} \begin{aligned} \overline {{\unicode{x1D654}}}\big (\boldsymbol{p}_2^R + \boldsymbol{p}_2^L\big ) &+ \overline {\mathcal{Y}}\big (\langle \boldsymbol{p}_2^R,\boldsymbol{p}_2^R\rangle + \langle \boldsymbol{p}_2^L,\boldsymbol{p}_2^L\rangle \big ) \\ &= {\unicode{x1D641}}^{{T}} \big ({\unicode{x1D654}}^{\,R}\boldsymbol{p}_1^R - {\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^L + \mathcal{Y}^R\langle \boldsymbol{p}_1^R,\boldsymbol{p}_1^R\rangle - \mathcal{Y}^L\langle \boldsymbol{p}_1^L,\boldsymbol{p}_1^L\rangle \big ), \end{aligned} \end{align}

from which $\boldsymbol{p}_2^L$ may be eliminated by imposing continuity of pressure in the outer duct away from the outlet; that is, by using the zero extension of the pressure discontinuity

(3.25) \begin{equation} \boldsymbol{p}_2^L = \boldsymbol{p}_2^R + {\unicode{x1D641}}^{{T}}\!\big (\boldsymbol{p}_1^L - \boldsymbol{p}_1^R\big ), \end{equation}

resulting in the equation for $\boldsymbol{p}_2^R$ in terms of $\boldsymbol{p}_1^R$ and $\boldsymbol{p}_1^L$

(3.26) \begin{align} \boldsymbol{p}_2^R = &\frac {1}{2}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}} \big ({\unicode{x1D654}}^{\,R}\boldsymbol{p}_1^R - {\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^L\big ) + \frac {1}{2}{\unicode{x1D641}}^{{T}} \big (\boldsymbol{p}_1^R - \boldsymbol{p}_1^L\big ) \notag \\ &\qquad + \frac {1}{2}\overline {{\unicode{x1D655}}}\Big \{{\unicode{x1D641}}^{{T}} \big (\mathcal{Y}^R\big \langle \boldsymbol{p}_1^R,\boldsymbol{p}_1^R\big \rangle - \mathcal{Y}^L\big \langle \boldsymbol{p}_1^L,\boldsymbol{p}_1^L\big \rangle \big ) \notag\\ &\qquad \qquad \qquad -\frac {1}{4}\overline {\mathcal{Y}}\Big (\big \langle \overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!\big ({\unicode{x1D654}}^{\,R}\boldsymbol{p}_1^R - {\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^L\big ) + {\unicode{x1D641}}^{{T}}\big (\boldsymbol{p}_1^R - \boldsymbol{p}_1^L\big ),\ldots \big \rangle \\ &\qquad \qquad \qquad \qquad \qquad + \big \langle \overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!\big ({\unicode{x1D654}}^{\,R}\boldsymbol{p}_1^R - {\unicode{x1D654}}^{\,L}\boldsymbol{p}_1^L\big ) - {\unicode{x1D641}}^{{T}}\big (\boldsymbol{p}_1^R - \boldsymbol{p}_1^L\big ),\ldots \big \rangle \Big )\Big \}. \notag \end{align}

Taking the restriction of both sides maps $\boldsymbol{p}_2^R$ to $\boldsymbol{p}_1^R$ ; we then rearrange for $\boldsymbol{p}_1^L$ , getting

(3.27) \begin{equation} \boldsymbol{p}_1^L = {\unicode{x1D63E}}\boldsymbol{p}_1^R + \big ({\unicode{x1D644}} + {\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,L}\big )^{-1}{\unicode{x1D641}}\mathcal{C} \big\langle \boldsymbol{p}_1^R,\boldsymbol{p}_1^R \big\rangle , \end{equation}

where $\unicode{x1D63E}$ and $\mathcal{C}$ are defined as

(3.28a) \begin{gather} {\unicode{x1D63E}} = -\big ({\unicode{x1D644}} + {\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,L}\big )^{-1} \big ({\unicode{x1D644}} - {\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,R}\big ), \end{gather}

(3.28b) \begin{gather} \mathcal{C} = \overline {{\unicode{x1D655}}}\left [{\unicode{x1D641}}^{{T}}\big (\mathcal{Y}^R - \mathcal{Y}^L\langle {\unicode{x1D63E}},{\unicode{x1D63E}}\rangle \big ) - \frac {1}{2}\overline {\mathcal{Y}}\Big (\big \langle \overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!\big ({\unicode{x1D654}}^{\,R} - {\unicode{x1D654}}^{\,L}{\unicode{x1D63E}}\big ),{\cdots} \big \rangle + \big \langle {\unicode{x1D641}}^{{T}}\!\big ({\unicode{x1D644}} - {\unicode{x1D63E}}\big ),{\cdots} \big \rangle \Big )\right ]\!. \end{gather}

Substituting this into the original equation (3.26) for $\boldsymbol{p}_2^R$ , we get

(3.29) \begin{align} \boldsymbol{p}_2^R = \frac {1}{2}\Big [\big (\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,R} + {\unicode{x1D641}}^{{T}}\big ) -& \big (\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,L} + {\unicode{x1D641}}^{{T}}\big ){\unicode{x1D63E}}\Big ]\boldsymbol{p}_1^R \nonumber \\ & + \frac {1}{2}\Big [{\unicode{x1D644}} - \big (\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,L} + {\unicode{x1D641}}^{{T}}\big )\big ({\unicode{x1D644}} + {\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,L}\big )^{-1}\!{\unicode{x1D641}}\Big ]\mathcal{C}\langle \boldsymbol{p}_1^R,\boldsymbol{p}_1^R\rangle , \end{align}

and then using (3.25) gives

(3.30) \begin{align} \boldsymbol{p}_2^L = \frac {1}{2}\Big [\big (\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,R} - {\unicode{x1D641}}^{{T}}\big ) -& \big (\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,L} - {\unicode{x1D641}}^{{T}}\big ){\unicode{x1D63E}}\Big ]\boldsymbol{p}_1^R \nonumber \\ &+ \frac {1}{2}\Big [{\unicode{x1D644}} - \big (\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,L} - {\unicode{x1D641}}^{{T}}\big )\big ({\unicode{x1D644}} + {\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}\!{\unicode{x1D654}}^{\,L}\big )^{-1}\!{\unicode{x1D641}}\Big ]\mathcal{C}\langle \boldsymbol{p}_1^R,\boldsymbol{p}_1^R\rangle . \end{align}

Equations (3.29) and (3.30) give the pressure at $s=0^+$ and $s=0^-$ as boundary conditions from which the pressure within the outer duct for $s\gt 0$ and $s\lt 0$ can be solved, respectively. The method of solution is the same as for the inner duct (again as described by Jensen & Brambley Reference Jensen and Brambley2026). This allows the acoustic field within the outer duct, modelling free space, to be plotted and visualised. Examples of this will be given below.

3.5. Numerical method

In this section we describe details of our numerical implementation of the open-end exit calculation outlined in § 3.3, and the subsequent pressure matching described in § 3.4. Matlab source code numerically implementing this model is included in the supplementary material; this code was used to produce the results that follow in § 5. For details of the numerical method for calculating propagation along either the inner or outer ducts, and a demonstration of its convergence as more modes are added, the reader is referred to Jensen & Brambley (Reference Jensen and Brambley2026).

When calculating the duct exit condition the infinite series of modes must be truncated during numerical implementation. Regarding truncation, it should be noted that the decision (taken in § 3.3) to set ${\unicode{x1D654}}_1^{\,L} = \overline {{\unicode{x1D654}}}_1$ and $\mathcal{Y}_1^L = \overline {\mathcal{Y}}_1$ on the acoustic wormhole has a numerical justification as well as a physical one. The natural alternative to characteristic admittances to the left of the outlet would be to take a zero admittance (corresponding to a hard wall and allowing no sound to propagate through): this causes numerical difficulties with the definitions of $\unicode{x1D64C}$ and $\widetilde {{\unicode{x1D64C}}}$ , and even when the equations are reformulated to avoid these problems, numerical difficulties still arise as a result of the non-invertibility of ${\unicode{x1D641}}\overline {{\unicode{x1D655}}}{\unicode{x1D641}}^{{T}}$ and ${\unicode{x1D641}}\overline {{\unicode{x1D654}}}{\unicode{x1D641}}^{{T}}$ when truncated. In contrast, for the perfectly absorbing characteristic admittance used here, no such numerical difficulties arise.

3.5.1. Choosing modenumbers and the width ratio

When truncating, we generally follow the rule of thumb set by McTavish (Reference McTavish2018) that $\alpha _{\textit{max}}^2 \geqslant \alpha _{\textit{max}}^1/\eta$ , where $\alpha _{\textit{max}}^1$ and $\alpha _{\textit{max}}^2$ are the maximum modenumbers in the inner and outer ducts, respectively. As well as choosing modenumbers carefully, it is also important that $\eta$ be small enough to approximate free space; there is therefore a trade-off between needing small $\eta$ , large $\alpha _{\textit{max}}^2$ and manageable runtime. The computations we perform with this framework are not always of the same complexity, since we will focus on different parts of the solution depending on what we wish to calculate. While the number of modes in the inner duct $\alpha _{\textit{max}}^1$ always remains around 10, $\alpha _{\textit{max}}^2$ and $\eta$ can vary much more.

Figure 2.

End-correction coefficients for a 3-D duct, calculated using the Wiener–Hopf technique (solid line) and our outlet condition (dash-dotted line), for $m_\alpha =0$ ; 8 modes were used in the inner duct and 1200 in the outer duct, with a duct width ratio of 1/40.

When only a single computation of the linear-regime exit condition at a given location is required (for example, when calculating the linear end correction described in § 4, or when calculating resonances in § 5.3), with no need to solve any ordinary differential equations (ODEs) in the outer duct, large numbers of modes can be used. In such cases, we will typically take of the order of a thousand modes for $\alpha _{\textit{max}}^2$ , meaning that we can take $\eta$ down to $1/40$ without under-resolving the duct exit. Such computations will only require tens of megabytes of memory in the linear case. These memory requirements apply in either two dimensions or three dimensions, but the calculation can take much longer in three dimensions than two dimensions, since the 3-D calculation involves integration of millions of Bessel function products: a 3-D calculation will often last hours on a standard desktop computer using a single core, whereas a 2-D one will be done in seconds. Figures 2 and 3 provide a means of calibrating $\eta$ and $\alpha _{\textit{max}}^2$ by validating against results from the Wiener–Hopf technique, as discussed further in § 5.1.3 below. We may still solve ODEs in the inner duct (i.e. in § 5.3) using this large value of $\alpha ^2_{\textit{max}}$ , since the number of modes in the inner duct $\alpha ^1_{\textit{max}}$ remains quite small, and linearity requires only a single temporal mode, so this does not present a time or memory constraint.

Figure 3.

End-correction coefficients for a 2-D duct with a symmetric source, calculated using the Wiener–Hopf technique (solid line) and our outlet condition (dash-dotted line). Thirteen modes were used in the inner duct and 2000 in the outer duct, with a width ratio of 1/40.

When nonlinearity becomes involved, memory requirements are a greater concern due to the presence of third-rank tensors. Using thousands of modes for a nonlinear calculation of the exit condition at a single location would require more than $128\,\mathrm{GB}$ of memory, so the number of modes was reduced in this case. It is demonstrated for the linear regime in § 5.1.3 that the parameters $(\eta ,\alpha _{\textit{max}}^1,\alpha _{\textit{max}}^2) = (1/10,10,200)$ still provide a reasonable approximation to the semi-analytical Wiener–Hopf solution: this requires of the order of $10\,\mathrm{MB}$ of memory even in the nonlinear case.

Computational time is also more of a limiting factor in the nonlinear regime (such as in §§ 5.2 and 5.6), especially since we will see in § 4.3 that the nonlinear exit condition cannot be considered in isolation without also solving for the acoustic field within the inner duct. In contrast to the linear regime, despite the number of spatial modes in the inner duct still being low, temporal modes are also required, and the ODEs are coupled in more ways than in the linear case. Performing a single computation for a simple duct geometry in the nonlinear case will thus take of the order of tens of minutes, meaning that the data for plots like figures 5 and 6 can take days to generate.

Moreover, even when calculating just the linear-acoustic field in the outer duct, we also need to reduce the number of modes relative to a single exit-condition calculation: using the number of modes from § 5.1 would cause calculations to run very slowly. In order to compute in the outer duct in a reasonable computation time, the choices of $(\alpha _{\textit{max}}^1,\alpha _{\textit{max}}^2,\eta ) = (10,200,1/10)$ are therefore maintained, resulting in calculations that take of the order of tens of seconds. This is further compounded in the nonlinear case: even with this more modest number of modes, combination with the $\sim 10$ temporal modes required for an accurate nonlinear calculation is beyond the memory and time constraints of our Matlab code on a standard desktop computer: we therefore restrict modelling of the nonlinear regime here to the inner-duct domain. While the emphasis in this paper is on the theoretical derivation and some proof-of-concept numerical examples to showcase possible behaviours, future work could focus on implementing the model in a more memory-efficient language than Matlab, or parallelising the code to run on a distributed-memory cluster, either of which would potentially make nonlinear outer-duct computations possible.

In summary, we use the numerical parameters $(\alpha ^1_{\textit{max}},\alpha _{\textit{max}}^2,\eta ) \!=\! (10,200,1/10)$ for the examples that follow, unless indicated otherwise. This is shown in § 5.1.3 to give sufficiently accurate results, although higher resolutions are used in a few special cases where they are needed for accuracy, as noted in the text.

4.1. Plane-wave case

If we first consider plane waves in the linear regime, the pressure within the inner duct is given by the superposition of an incident plane wave propagating towards the exit at $s_{\text{o}}$ and a reflected plane wave propagating back towards the entrance with complex reflection coefficient $R$

(4.1) \begin{equation} p = \text{Re}\Big [p_0\big ({\mathrm{e}}^{\gamma s - {i}\omega t} - |R|{\mathrm{e}}^{-\gamma s - {i}\omega t + {i}\arg (-R)}\big )\Big ], \end{equation}

with $-\mathrm{i}\gamma \gt 0$ (since plane waves are always cut-on). Assuming that a node in the pressure occurs at some distance beyond the outlet $\hat {s}$ $ = s - s_{\text{o}}$ , and taking the node to be a point where the incident and reflected waves are in phase and so maximally cancel, we get

(4.2) \begin{align} &\mathrm{Im}\big [\gamma \hat {s}-\mathrm{i}\omega t\big ] = \mathrm{Im}\big [-\gamma \hat {s}-\mathrm{i}\omega t + \mathrm{i}\arg (-R)\big ] & &\Rightarrow & \hat {s} &= \frac {\arg (-R)}{2\mathrm{Im}(\gamma )}. \end{align}

4.2. Linear generalisation

Attempting to generalise the above for non-plane-wave duct modes, we express the pressure as an expansion in the straight-duct modes at a single frequency (effectively setting $a=1$ in the language of § 2)

(4.3) \begin{equation} p = \text{Re}\left [{\mathrm{e}}^{-{i}\omega t}\sum _{\alpha }\psi _\alpha (\boldsymbol{x})\big ({\mathrm{e}}^{\gamma _\alpha s}A_{\alpha }^{+} + {\mathrm{e}}^{-\gamma _\alpha s}A_{\alpha }^{-}\big )\right]\!, \end{equation}

where $\gamma _\alpha$ is the straight-duct eigenvalue $\overline {\gamma }_\alpha$ from Jensen & Brambley (Reference Jensen and Brambley2026), given e.g. in three dimensions by

(4.4) \begin{equation} \gamma _\alpha = \begin{cases} \mathrm{i}\omega \sqrt {\big |1 - \lambda _\alpha ^2/\omega ^2R^2\big |}\quad &\text{Re}(\omega ) \geqslant \lambda _\alpha /R,\\[11pt] -\omega \sqrt {\big |1 - \lambda _\alpha ^2/\omega ^2R^2\big |} &\text{Re}(\omega ) \lt \lambda _\alpha /R, \end{cases} \end{equation}

and $A_\alpha ^+$ and $A_\alpha ^-$ are vectors of constants. The reflection matrix is given by

(4.5) \begin{equation} {\unicode{x1D64D}} = \big ({\unicode{x1D654}} + \overline {{\unicode{x1D654}}}\big )^{\!-1} \big (\overline {{\unicode{x1D654}}} - {\unicode{x1D654}}\big ), \end{equation}

and maps from the forward-going pressure to the backward, i.e. $\boldsymbol{p}^- = {\unicode{x1D64D}}\boldsymbol{p}^+$ . From the definition of the reflection matrix at the duct outlet, we may write the reflected waves in terms of the outgoing waves

(4.6) \begin{equation} A_\alpha ^{-}{\mathrm{e}}^{-\gamma _\alpha s_{\text{o}}} = P_\alpha ^{-}(s_{\text{o}}) = \sum _\beta {\unicode{x1D64D}}_{\alpha \beta }(s_{\text{o}})A_\beta ^{+}{\mathrm{e}}^{\gamma _\beta s_{\text{o}}}. \end{equation}

Writing $\hat {s} = s - s_{\text{o}}$ , and factoring out the $A_\beta ^{+}$ vector, we have an expression equivalent to (4.1)

(4.7) \begin{equation} p = \text{Re}\left[{\mathrm{e}}^{-{i}\omega t}\sum _{\alpha ,\beta } \psi _\alpha (\boldsymbol{x})\big ({\mathrm{e}}^{\gamma _\alpha \hat {s}}\delta _{\alpha \beta } - |{\unicode{x1D64D}}_{\alpha \beta }(s_{\text{o}})|{\mathrm{e}}^{-\gamma _\alpha \hat {s} + {i}\arg (-{\unicode{x1D64D}}_{\alpha \beta }(s_{\text{o}}))}\big )A_\beta ^{+}{\mathrm{e}}^{\gamma _\beta s_{\text{o}}}\right]\!; \end{equation}

since we do not have a sensible definition of the end correction for cut-off waves, we restrict ourselves to the cut-on case. Further to this, because we are comparing the arguments of entries in the reflection matrix with those of the matrix ${\mathrm{e}}^{2\gamma _\alpha \hat {s}}\delta _{\alpha \beta }$ , which is diagonal, we correspondingly only consider the reflection matrix’s diagonal entries. This results in a different end correction $\hat {s}_{\alpha }$ for each mode $\alpha$ that is cut-on (with $\mathrm{Im}(\gamma _\alpha )\neq 0$ ), given by

(4.8) \begin{equation} \hat {s}_{\alpha } = \frac {\arg (-{\unicode{x1D64D}}_{\alpha \alpha }(s_{\text{o}}))}{2\mathrm{Im}(\gamma _\alpha )}. \end{equation}

Given that we do not expect the oscillations in a musical instrument to be monomodal, but rather to consist of a series of harmonics, a set of modal end corrections will be informative as to how each harmonic’s pitch will be affected according to the reflective properties of the instrument’s outlet. This in turn will affect the timbre of the noise produced.

4.3. Nonlinear generalisation

It is not obvious how best to introduce nonlinearity to the end-correction problem. The reflection matrix $\unicode{x1D64D}$ may be generalised to weak nonlinearity exactly as for the admittance or impedance, giving $\boldsymbol{p}^- = {\unicode{x1D64D}}\boldsymbol{p}^+ + \mathcal{R}\langle \boldsymbol{p}^+,\boldsymbol{p}^+\rangle$ , but this reflects incoming modes into different reflected modes, and so does not directly affect the end correction; instead, these other reflected modes feed back on the initial mode through weakly nonlinear propagation along the inner duct. Here, we consider a different nonlinear correction to the reflection matrix that includes both nonlinear scattering at the open end and nonlinear interactions within the duct; this alternative nonlinear reflection matrix then facilitates the calculation of a nonlinear end correction. The linear reflection matrix at the open end, as given above in terms of the admittance, may also be obtained directly from the acoustic pressure in the duct if a single duct mode $(a,\alpha )$ is excited at the inlet, so that $P_\beta ^{b+} = A_\alpha ^{a+}\delta _{\alpha \beta }\delta ^{ab}$ . We take $a = 1$ in order to see how a source at the fundamental excites higher harmonics. In this case, at $s = s_{\mathrm{o}}$ , we have

(4.9) \begin{equation} \frac {\sum _{b = 1}^\infty \sum _\beta \big(P_\beta ^{b+}\big)^*P_\beta ^{b-}}{\sum _{b = 1}^\infty \sum _\beta |P_\beta ^{b+}|^2} = \frac {\sum _{b = 1}^\infty \sum _{\beta ,\gamma } \big(P_\beta ^{1+}\big)^*{\unicode{x1D64D}}_{\beta \gamma }^{\,b}P_\gamma ^{1+}}{\sum _{b = 1}^\infty \sum _\beta |P_\beta ^{1+}|^2} +O(M) = {\unicode{x1D64D}}^1_{\alpha \alpha } + O(M), \end{equation}

where ‘ $*$ ’ denotes a complex conjugate and the nonlinear reflection tensor $\mathcal{R}$ is contained within the $O(M)$ term. We may therefore define a matricial reflection coefficient $\tilde {{\unicode{x1D64D}}}$ that includes both the nonlinear propagation in the duct and the nonlinear reflection at the open end by taking

(4.10) \begin{align} & \tilde {{\unicode{x1D64D}}}_{\alpha \alpha }^1 := \frac {\sum _{b = 1}^\infty \sum _\beta \big(P_\beta ^{b+}\big)^*P_\beta ^{b-}}{\sum _{b = 1}^\infty \sum _\beta |P_\beta ^{b+}|^2}, &&\Rightarrow && \tilde {s}_{\alpha } = \frac {\arg \big(-\tilde {{\unicode{x1D64D}}}_{\alpha \alpha }^1\big)}{2\mathrm{Im}(\gamma _\alpha ^1)}, \end{align}

where $\tilde {s}_\alpha$ is a modified, nonlinear end correction. Note that this definition only holds if the fundamental $a = 1$ is the only mode excited at the inlet.

An important difference between this calculation and the linear one is that the length and geometry of the duct now matter (since these affect the pressure at the outlet). The amount of steepening that a wave can undergo within the duct will influence its behaviour at the outlet, so the previously 1-D parameter space (frequency only) is now three-dimensional (frequency, Mach number and duct length). We choose to fix the frequency and study the 2-D relationship between Mach number and duct length for the fundamental frequency ( $a=1$ ) plane-wave ( $\alpha = 0$ ) coefficient $\tilde {s}_0$ ; the two other combinations and higher mode entries could be considered in future work. Results from this calculation will be shown in § 5.2.