2.1. Constructing the Pridmore-Brown equation

The governing equations for what follows are the Euler equations with a mass source $q$:

(2.1a–c)\begin{equation} \frac{\partial {\rho} }{\partial {t} }+\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho \boldsymbol{u})=q, \quad \rho\,\frac{{\rm D} {\boldsymbol{u}} }{{\rm D} {t} }={-}\boldsymbol{\nabla}p, \quad \frac{{\rm D} {p} }{{\rm D} {t} }={c^{2}}\,\frac{{\rm D} {\rho} }{{\rm D} {t} }, \end{equation}

where $\boldsymbol {u}$ is the fluid velocity, $p$ is the pressure, $\rho$ is the density, and $c^{2} = {\partial p}/{\partial \rho }|_s$ is the square of the sound speed. In what follows, we take the mass source $q$ to be a small time-harmonic point mass source. In cylindrical coordinates $(x,r,\theta )$, with a suitable choice of origin, this mass source $q$ may in general be taken as

(2.2)\begin{equation} q = {\rm Re}\left(\frac{\epsilon}{r_0}\,\delta(x)\,\delta(\theta)\, \delta(r-r_0)\exp\{\mathrm{i}\omega t\}\right), \end{equation}

where $\epsilon$ is the small amplitude, $\omega$ is the frequency, and the $1/r_0$ term comes from writing a unit amplitude point source in cylindrical coordinates. We expand each variable in powers of $\epsilon$:

(2.3)\begin{equation} \left.\begin{array}{c@{}} \rho=\rho_0(r)+{\rm Re}(\epsilon\hat{\rho}\,{\rm e}^{\mathrm{i}\omega t}) + O(\epsilon^{2}), \quad p=p_0+{\rm Re}(\epsilon \hat{p}\, {\rm e}^{\mathrm{i}\omega t}) + O(\epsilon^{2}),\\ \boldsymbol{u}=U(r)\,\boldsymbol{e_x} + {\rm Re}(\epsilon (\hat{u}, \hat{v}, \hat{w} )\, {\rm e}^{\mathrm{i}\omega t}) + O(\epsilon^{2}), \quad c^{2} = c_0^{2}(r) + O(\epsilon), \end{array}\right\} \end{equation}

where $p_0$ is necessarily a constant in order that the steady state should satisfy the Euler equations, and it turns out that $c^{2}$ is needed only to leading order in what follows. Without loss of generality, all perturbations are expanded using a Fourier series in $\theta$ and a Fourier transform in $x$. As a result, the pressure perturbation is given as

(2.4)\begin{equation} \hat{p}(x,r,\theta)=\frac{1}{2{\rm \pi}} \sum_{m={-}\infty}^{\infty} \,\mathrm{e}^{-\mathrm{i} m \theta} \int_{-\infty}^{\infty} \tilde{p}(r;k,m,\omega)\, \mathrm{e}^{-\mathrm{i} k x}\,\mathrm{d} k, \end{equation}

and similarly for the density $\hat {\rho }$ and the velocity components $\hat {u}$, $\hat {v}$ and $\hat {w}$. Substituting these into the Euler equations (2.1a–c), and linearizing by ignoring terms of $O(\epsilon ^{2})$ or smaller, each of $\tilde {\rho }$, $\tilde {u}$, $\tilde {w}$ and finally $\tilde {v}$ may be eliminated, to leave a second-order ordinary differential equation in the radial coordinate $r$ for $\tilde {p}$:

(2.5a)\begin{gather} \tilde{p}^{\prime\prime}+\left( {\frac{2k U^{\prime}}{\omega-U(r)\,k}+ \frac{1}{r}-\frac{\rho^{\prime}_0}{\rho_0}} \right)\tilde{p}^{\prime} + \left( {\frac{(\omega-U(r)\,k)^{2}}{c_0^{2}}-k^{2}-\frac{m^{2}}{r^{2}}} \right) \tilde{p}\nonumber\\ \quad =\frac{\omega-U(r_0)\,k}{2\mathrm{i}{\rm \pi} r_0}\,\delta(r-r_0), \end{gather}

(2.5b)\begin{gather}\text{with}\quad\tilde{v} = \frac{\mathrm{i}\tilde{p}'}{\rho_0(\omega - Uk)}, \end{gather}

where a prime denotes the derivative with respect to $r$. This is the Pridmore-Brown (Reference Pridmore-Brown1958) equation for a point mass source, written in cylindrical coordinates.

One boundary condition to (2.5) is regularity at $r=0$. The singular solution behaves, for $m\neq 0$, as $O(r^{-|m|})$ as $r\to 0$, and the regular solution behaves as $O(r^{|m|})$. For $m=0$, the singular solution behaves as $O(\log r)$, while the regular solution behaves as $O(1)$. Eliminating the singular solution is therefore possible using the boundary conditions at $r=0$:

(2.6a,b)\begin{equation} \tilde{p}(0) = 0 \text{ for }m\not=0, \quad \tilde{p}^{\prime}(0)=0 \text{ for }m=0. \end{equation}

To model sound within a straight cylindrical duct of radius $r=a$, we take the other boundary condition to be the impedance boundary condition at $r=a$:

(2.7)\begin{equation} \tilde{p}({a}) = Z(\omega)\,\tilde{v}({a}) \iff \tilde{p}^{\prime}({a})={-}\frac{\mathrm{i}\omega{\rho_0}}{Z}\,\tilde{p}({a}), \end{equation}

where $Z(\omega )$ is the impedance of the duct wall, and the two expressions are equivalent in light of (2.5b). A hard wall corresponds to $Z \to \infty$, and hence to $\tilde {v}(a) = 0$, or equivalently to $\tilde {p}'(a) = 0$.

In what follows, we make the simplifying assumption of a constant density $\rho _0(r)$. This is a homentropic assumption, and it implies that $c_0(r)$ is also constant. We may then non-dimensionalize speeds by the sound speed $c_0$, densities by $\rho _0$, and distances by the duct radius $a$. Note that this places the impedance boundary condition in non-dimensional terms at $r=1$. We also assume a flow profile $U(r)$ that is uniform, except within a boundary layer of width $h$ where it varies quadratically:

(2.8)\begin{equation} U(r)=\begin{cases}M, & 0\leq r\leq 1-h, \\ M\left(1-\left(1-\dfrac{1-r}{h}\right)^{2}\right), & 1-h\leq r \leq 1.\end{cases} \end{equation}

With the non-dimensionalization of velocities by $c_0$, $M$ here is the duct centreline Mach number. This situation is depicted schematically in figure 1.

In order to solve the Pridmore-Brown equation (2.5a), we first consider solutions of the homogeneous form

(2.9)\begin{equation} \tilde{p}^{\prime\prime}+\left( {\frac{2k U^{\prime}}{\omega-U(r)\,k}+ \frac{1}{r}} \right)\tilde{p}^{\prime}+\left( {(\omega-U(r)\,k)^{2}-k^{2}-\frac{m^{2}}{r^{2}}} \right)\tilde{p}=0. \end{equation}

2.2. Homogeneous solutions within the region of uniform flow

Within the region of uniform flow, the homogeneous Pridmore-Brown equation (2.9) reduces to

(2.10)\begin{equation} \tilde{p}^{\prime\prime}+\frac{1}{r}\tilde{p}^{\prime}+\left( {(\omega-Mk)^{2}-k^{2}- \frac{m^{2}}{r^{2}}} \right)\tilde{p}=0. \end{equation}

This is Bessel's equation of order $m$ rescaled by $\alpha$, where

(2.11)\begin{equation} \alpha^{2}=(\omega-Mk)^{2}-k^{2}; \end{equation}

it will turn out later that the branch chosen for $\alpha$ does not matter, although for definiteness, one may choose $\textrm {Re}(\alpha )>0$. Bessel's equation has two pairs of linearly independent solutions that we will make use of: the Bessel functions of the first and second kind, $J_m(\alpha r)$ and $Y_m(\alpha r)$; and the Hankel functions of the first and second kind, $H_m^{(1)}(\alpha r)$ and $H_m^{(2)}(\alpha r)$. More information regarding these can be found in Abramowitz & Stegun (Reference Abramowitz and Stegun1964). It is worth noting that only $J_m(\alpha r)$ is regular at $r=0$, with the other solutions all requiring a branch cut along $\alpha r<0$, with a singularity at $\alpha r=0$.

2.3. Homogeneous solutions within the region of sheared flow

In this section, we will construct the solution of the homogeneous Pridmore-Brown equation (2.9) when $U(r)$ varies by proposing a Frobenius expansion about the singularities of the Pridmore-Brown equation.

In addition to the singularity at $r=0$, the homogeneous Pridmore-Brown equation possesses regular singularities whenever $\omega - U(r)\,k = 0$; these singularities correspond to the critical layer. Within the sheared flow region $1-h< r<1$, since the velocity profile $U(r)$ is quadratic in $r$, there are exactly two critical values $r = r_c$ for which $\omega -U(r_c)\,k=0$. Note that in general, these critical values will be complex. Solving this quadratic equation gives the two singularities explicitly as $r_c^{+}$ and $r_c^{-}$, where

(2.12)\begin{equation} r_c^{{\pm}}=1-h\pm Q, \quad Q=h\sqrt{1-\frac{\omega}{Mk}}. \end{equation}

For convenience, we will take $\textrm {Re}(Q)\geq 0$, so that $\textrm {Re}(r_c^{+}) \geq 1-h$ and $\textrm {Re}(r_c^{-}) \leq 1-h$. Since solutions with this quadratic flow profile $U(r)$ are valid only for $1-h< r<1$, it will therefore be $r_c^{+}$ that we are mostly concerned about here.

Following Brambley et al. (Reference Brambley, Darau and Rienstra2012a), we propose a Frobenius expansion (Teschl Reference Teschl2012) about the regular singularity $r_c^{+}$:

(2.13)\begin{equation} \tilde{p}(r)=\sum_{n=0}^{\infty} a_n(r-r_c^{+})^{n+\sigma}, \quad \text{with } a_0 \neq 0. \end{equation}

Specifying $a_0 \neq 0$ results in a condition on $\sigma$, and we find that $\sigma =0,3$. By Fuchs’ theorem (Teschl Reference Teschl2012), this gives a pair of linearly independent solutions of the form

(2.14a)\begin{gather} \tilde{p}_1(r)=\sum_{n=0}^{\infty} a_n (r-r_c^{+})^{n+3}, \end{gather}

(2.14b)\begin{gather}\tilde{p}_2(r)=A\,\tilde{p}_1(r)\log(r-r_c^{+})+\sum_{n=0}^{\infty} b_n(r-r_c^{+})^{n}. \end{gather}

The coefficients $a_n$ and $b_n$ are derived in Appendix A, where, in particular, it is found that

(2.15a–d)\begin{equation} a_0 = b_0 = 1, \quad b_1 = 0, \quad b_2 ={-}\frac{1}{2} \left(k^{2}+\left(\frac{m}{r_c^+}\right)^2\right) , \quad b_3 =0, \end{equation}

and that

(2.16)\begin{equation} A={-}\frac{1}{3}\left( {\frac{1}{Q}-\frac{1}{r_c^{+}}} \right) \left(k^{2} + \left(\frac{m}{r_c^+}\right)^2\right) -\frac{2m^{2}}{3r_c^{{+}3}}, \end{equation}

the latter in agreement with (2.3)–(2.5) of Brambley et al. (Reference Brambley, Darau and Rienstra2012a). We note in passing that in practice we may be limited by the radius of convergence of (2.14), and in such cases, the solutions given above are analytically continued by a companion expansion of the Pridmore-Brown equation about $r=1$, as described in § A.2. Other than being a complication concerning numerical convergence, this complication may be ignored, and $\tilde {p}_1$ and $\tilde {p}_2$ thought of as being defined by the expressions in (2.14).

Due to the log term in $\tilde {p}_2$ in (2.14b), a branch cut is necessary in the complex $r$-plane originating from the branch point $r = r_c^{+}$. This branch cut must be such that the solutions remain continuous for the real values of $r\in [1-h,1]$, so the branch cut must avoid crossing the real $r$-axis between $1-h$ and $1$. In the following, we achieve this by choosing the branch cuts parallel to the imaginary axis and away from the real axis, as depicted in figure 2. When $r_c^{+}$ is real and $1-h< r_c^{+}<1$, no suitable choice of branch cut exists, and as a result, any solution $\tilde {p}(r)$ with $\tilde {p}(r_c^{+})\neq 0$ necessarily has a singular third derivative at $r_c^{+}$. This occurs only for particular values of $k$, however, and we can map the corresponding values of $k$ in the complex $k$-plane to find that they fall exactly on the half-line $[{\omega }/{M},\infty )$; we refer to this range of excluded values of $k$ as the critical layer branch cut. As $r_c^{+}$ becomes real, note that the value of $\tilde {p}_2(r_c^{+})$ is different depending on whether we approach from the positive or negative imaginary part. Thinking of $r_c^{+}(k)$ as a function of $k$, this corresponds to approaching the critical layer branch cut $[{{\omega }/{M}},\infty )$ in $k$ from above or below. This reinforces the consideration of the critical layer appearing as a branch cut in the complex $k$-plane. The change in $\tilde {p}_2$ when crossing the critical layer branch cut from below to above is described as

(2.17)\begin{equation} {\rm \Delta} \tilde{p}_2(r) = \lim_{{{{\rm Im}}(k)\searrow 0+}} \tilde{p}_2(r) - \lim_{{{{\rm Im}}(k)\nearrow 0-}} \tilde{p}_2(r) ={-}2{\rm \pi} \mathrm{i} A\,\tilde{p}_1(r)\,H(r_c^{+}-r), \end{equation}

where $H(r)$ is the Heaviside function.

Figure 2.

Schematic of possible locations of the $r_c^{+}$ branch cut in the complex $r$-plane. (a) A possible choice of branch cut when $\textrm {Im} (r_c^{+})>0$. (b) The other choice of branch cut is needed when $\textrm {Im} (r_c^{+})<0$.

In order to retrieve this result, we need only consider the $\log (r-r_c^{+})$ term of $\tilde {p}_2$. Note that $\partial r_c^{+}/\partial k > 0$ for real $k$ and real positive $\omega$; hence if $k$ is nearly real and $\textrm {Im} (k)>0$, then $\textrm {Im} (r_c^{+})>0$, and we must take the branch cut of $\log (r-r_c^{+})$ upwards towards $+\mathrm {i}\infty$. Similarly, if $\textrm {Im} (k)<0$, then $\textrm {Im} (r_c^{+})<0$, and the branch cut for $\log (r-r_c^{+})$ must be taken downwards to $-\mathrm {i}\infty$. When $r_c^{+}$ is located on the real line, $(r-r_c^{+})$ is negative for $r< r_c^{+}$. When we choose the branch cut into the upper half-plane, this corresponds to a complex argument of $-{\rm \pi}$. When we choose the branch cut into the lower half-plane, this corresponds to a complex argument of ${\rm \pi}$. This difference results in the jump of $2{\rm \pi} \mathrm {i}$ given. If we instead consider $r>r_c^{+}$, then the same argument is retrieved regardless of which direction we take the branch cuts, so no jump is observed. This is the reason for the presence of the Heaviside function.

2.4. Homogeneous solutions across the full domain

In order to construct a full solution in $r\in [0,1]$, we now construct two solutions $\psi _1(r)$ and $\psi _2(r)$ that solve (2.9) across $r\in [0,1]$, by patching together the solutions derived above in §§ 2.2 and 2.3. We construct $\psi _1(r)$ to satisfy the boundary conditions (2.6a,b) at $r=0$, by taking

(2.18)\begin{equation} \psi_1(r)=\begin{cases} J_m(\alpha r), & 0\leq r\leq 1-h, \\ C_1\,\tilde{p}_1(r)+D_1\,\tilde{p}_2(r), & 1-h\leq r\leq 1, \end{cases} \end{equation}

where the coefficients $C_1$ and $D_1$ ensure $C^{1}$ continuity, and are given by

(2.19a)\begin{gather} C_1=\frac{J_m(\alpha(1-h))\,\tilde{p}_2^{\prime}(1-h)- \alpha\,J_m^{\prime}(\alpha(1-h))\,\tilde{p}_2(1-h)}{W(1-h)}, \end{gather}

(2.19b)\begin{gather}D_1={-}\frac{J_m(\alpha(1-h))\,\tilde{p}_1^{\prime}(1-h)-\alpha\, J_m^{\prime}(\alpha(1-h))\,\tilde{p}_1(1-h)}{W(1-h)}, \end{gather}

and $W(r) = \mathcal {W}(\tilde {p}_1,\tilde {p}_2;r)$ is the Wronskian of $\tilde {p}_1$ and $\tilde {p}_2$, given in § A.4 as

(2.20)\begin{equation} W(r)=\mathcal{W}(\tilde{p}_1,\tilde{p}_2;r) = \tilde{p}_1(r)\,\tilde{p}_2'(r)-\tilde{p}_{2}(r)\,\tilde{p}_1'(r) ={-}\frac{3}{4}\,\frac{r_c^{+}(r-r_c^{+})^{2}(r-r_c^{-})^{2}}{rQ^{2}}. \end{equation}

Having constructed $\psi _1$ to satisfy the boundary condition at $r=0$, we now proceed to construct $\psi _2$ that satisfies the boundary condition (2.7) at $r=1$. Writing $\psi _2$ in terms of the homogeneous solutions derived above,

(2.21)\begin{equation} \psi_2(r)=\begin{cases} \check{C}_2\,H_m^{(1)}(\alpha r)+\check{D}_2\,H_m^{(2)}(\alpha r), & 0\leq r\leq 1-h, \\ \hat{C}_2\,\tilde{p}_1(r)+\hat{D}_2\,\tilde{p}_2(r), & 1-h\leq r \leq 1, \end{cases} \end{equation}

we choose $\hat {C}_2$ and $\hat {D}_2$ to satisfy $\psi _2(1)=1$ and $\psi ^{\prime }_2(1)=-{\mathrm {i}\omega }/{Z}$. This forces a non-zero normalized solution for $\psi _2$ that satisfies the boundary condition (2.7) at $r=1$, and leads to

(2.22a,b)\begin{equation} \hat{C}_2= \frac{{\tilde{p}_2'(1)}+\dfrac{\mathrm{i}\omega}{Z}\,{\tilde{p}_2(1)}}{W(1)}, \quad \hat{D}_2={-}\frac{{\tilde{p}_1'(1)}+\dfrac{\mathrm{i}\omega}{Z}\,{\tilde{p}_1(1)}}{W(1)}. \end{equation}

The coefficients $\check {C}_2$ and $\check {D}_2$ are chosen such that our solution is $C^{1}$ continuous at $r=1-h$, giving

(2.23)\begin{align} \left( { \begin{array}{c} \check{C}_2 \\ \check{D}_2 \end{array} } \right) &=\frac{\mathrm{i}{\rm \pi}(1 - h)}{4} \left( { \begin{array}{cc} \alpha\,H_m^{(2)\prime} (\alpha (1 - h)) & {-}H_m^{(2)}(\alpha (1 - h)) \\ {-}\alpha\,H_m^{(1)\prime}(\alpha (1 - h)) & H_m^{(1)}(\alpha (1 - h)) \end{array} } \right)\nonumber\\ &\quad\times \left( { \begin{array}{cc} \tilde{p}_1(1 - h) & \tilde{p}_2(1 - h) \\ \tilde{p}_1^{\prime}(1 - h) & \tilde{p}_2^{\prime}(1 - h) \end{array} } \right) \left( { \begin{array}{c} \hat{C}_2 \\ \hat{D}_2 \end{array} } \right), \end{align}

where the factor at the beginning comes from the Wronskian of $H_m^{(1)}$ and $H_m^{(2)}$ from Abramowitz & Stegun (Reference Abramowitz and Stegun1964), formula 9.1.17.

We will also require later the jump in behaviour of $\psi _1$ and $\psi _2$ as $k$ crosses the critical layer branch cut from below to above. Since any jump comes from the log term in $\tilde {p}_2(r)$ when $r< r_c^{+}$, we have, provided that $r_c^{+}<1$,

(2.24a)\begin{gather} {\rm \Delta} C_1=2\mathrm{i}{\rm \pi} A D_1, \quad {\rm \Delta} \hat{C}_2 = {\rm \Delta} D_1 = {\rm \Delta} \hat{D}_2 = 0, \end{gather}

(2.24b)\begin{gather}\left( { \begin{array}{c} {\rm \Delta} \check{C}_2 \\ {\rm \Delta} \check{D}_2 \end{array} } \right)= \frac{{\rm \pi}^{2}(1 - h)A \hat{D}_2}{2} \left( { \begin{array}{cc} \alpha\,H_m^{(2)\prime} (\alpha (1 - h)) & {-}H_m^{(2)}(\alpha (1 - h)) \\ {-}\alpha\,H_m^{(1)\prime}(\alpha (1 - h)) & H_m^{(1)}(\alpha (1 - h)) \end{array}} \right) \left( { \begin{array}{c} \tilde{p}_1(1 - h) \\ \tilde{p}_1^{\prime}(1 - h) \end{array}} \right), \end{gather}

resulting in (provided that $r_c^{+}<1$)

(2.25a)\begin{gather} {\rm \Delta} \psi_1(r)=2\mathrm{i}{\rm \pi} AD_1\,\tilde{p}_1H(r-r_c^{+}), \end{gather}

(2.25b)\begin{gather}{\rm \Delta} \psi_2(r)=\begin{cases} {\rm \Delta} \check{C}_2\,H_m^{(1)}(\alpha r)+{\rm \Delta} \check{D}_2\,H_m^{(2)}(\alpha r), & 0\leq r\leq 1-h, \\ -2\mathrm{i}{\rm \pi} A\,\tilde{p}_1(r)\,\hat{D}_2\,H(r_c^{+}-r), & 1-h\leq r \leq 1. \end{cases} \end{gather}

Note that if $r_c^{+}>1$, then ${\rm \Delta} \psi _1 = {\rm \Delta} \psi _2 = 0$, since the $\psi _1$ and $\psi _2$ solutions are defined uniquely by their boundary conditions, and no branch point occurs on the interval $r\in [1-h,1]$ to cause a jump.

2.5. Modal solutions

Modal solutions of the homogeneous Pridmore-Brown equation (2.9) are non-zero solutions $\tilde {p}(r)$ that satisfy the boundary conditions (2.6a,b) and (2.7) at $r=0$ and $r=1$. In general, satisfying both boundary conditions would force the solution $\tilde {p}(r) \equiv 0$, so non-zero solutions exist only for particular modal eigenvalues $k$ (assuming that $\omega$ is given and fixed). In contrast, the solution $\psi _1(r)$ is never identically zero and always satisfies the homogeneous Pridmore-Brown equation and the boundary condition at $r=0$; indeed, any solution satisfying the boundary condition at $r=0$ is necessarily a multiple of $\psi _1(r)$. Likewise, the solution $\psi _2(r)$ is never identically zero and always satisfies the homogeneous Pridmore-Brown equation and the boundary condition at $r=1$, and any solution satisfying the boundary condition at $r=1$ is necessarily a multiple of $\psi _2(r)$. In general, $\psi _1$ and $\psi _2$ are linearly independent, so their Wronskian $\mathcal {W}(\psi _1,\psi _2;r)$ is not identically zero, where

(2.26)\begin{equation} \mathcal{W}(\psi_1,\psi_2;r)=\psi_1(r)\,\psi_2^{\prime}(r)- \psi_2(r)\,\psi_1^{\prime}(r). \end{equation}

However, if $\tilde {p}(r)$ is non-zero and satisfies both the boundary conditions at $r=0$ and $r=1$, then $\tilde {p}(r) = a\,\psi _1(r) = b\,\psi _2(r)$ for some non-zero coefficients $a, b$. In other words, a modal solution is one where $\psi _1$ and $\psi _2$ are linearly dependent, so $\mathcal {W}(\psi _1,\psi _2;r) \equiv 0$.

For $1-h\leq r\leq 1$, substituting $\psi _1$ from (2.18) and $\psi _2$ from (2.21) into the Wronskian (2.26) gives

(2.27)\begin{equation} \mathcal{W}(\psi_1,\psi_2;r)=(C_1\hat{D}_2-\hat{C}_2D_1)\,W(r), \end{equation}

where $W(r)$ is the Wronskian between $\tilde {p}_1$ and $\tilde {p}_2$ and is given earlier, in (2.20). Since $\tilde {p}_1$ and $\tilde {p}_2$ were constructed to be linearly independent, we expect $W(r)$ not to be identically zero, and indeed (2.20) shows that $W(r) \neq 0$ except at the critical layer $r = r_c^{+}$. A modal solution, therefore, is given by the condition $C_1\hat {D}_2 - \hat {C}_2 D_1 = 0$, which is independent of $r$, and implies that $C_1/D_1 = \hat {C}_2/\hat {D}_2$, so that $\psi _1$ and $\psi _2$ are multiples of one another.

The same can be seen for $r\leq 1-h$. In this case, the Wronskian (2.26) becomes

(2.28)\begin{equation} \mathcal{W}(\psi_1,\psi_2;r) = \alpha\check{C}_2\,\mathcal{W}(J_m,H_m^{(1)};r) +\alpha\check{D}_2\,\mathcal{W}(J_m,H_m^{(2)};r) = \alpha(\check{C}_2-\check{D}_2)\,\frac{2\mathrm{i}}{{\rm \pi} r}, \end{equation}

where we have made use of the Bessel function identities 9.1.3, 9.1.4 and 9.1.16 from Abramowitz & Stegun (Reference Abramowitz and Stegun1964). Note in particular that $r\,\mathcal {W}(\psi _1,\psi _2;r)$ is a constant independent of $r$ for $0\leq r\leq 1-h$. Since $\mathcal {W}(\psi _1,\psi _2;r)$ is continuous in $r$ across $r = 1-h$, because $\psi _1$ and $\psi _2$ are both $C^{1}$ continuous, it follows that for $0\leq r\leq 1-h$, we can set $r\,\mathcal {W}(\psi _1,\psi _2;r)=(1-h)\,\mathcal {W}(\psi _1,\psi _2;1-h)$. We therefore arrive at the conclusion that

(2.29)\begin{equation} \mathcal{W}(\psi_1,\psi_2;r) = (C_1\hat{D}_2-\hat{C}_2D_1)\times\begin{cases} W(r), & 1-h \leq r \leq 1, \\ W(1-h)\,\dfrac{1-h}{r}, & 0 \leq r \leq 1-h, \end{cases} \end{equation}

and that a mode corresponds to the dispersion relation $0 = D(k,\omega ) = C_1\hat {D}_2-\hat {C}_2D_1$. In the next section, we see how these modal solutions occur naturally as poles in the solution of the non-homogeneous Pridmore-Brown equation.

3.1. Inhomogeneous solution of the Pridmore-Brown equation

While previously we have been solving only the homogeneous form (2.9), our original problem was to solve the inhomogeneous Pridmore-Brown equation (2.5a) subject to a harmonic point mass source. Due to the right-hand side of (2.5a) being a scalar multiple of a delta function, located at $r=r_0$, our solution will be the same scalar multiple of the Green's function, and we denote this solution as $\tilde {G}$. This function will satisfy the boundary condition at $r=0$ and $r=1$, and will solve the homogeneous Pridmore-Brown equation for $r< r_0$ and $r>r_0$; hence $\tilde {G}$ may be written as a multiple of the homogeneous solution $\psi _1$ for $r< r_0$, and as a multiple of the homogeneous solution $\psi _2$ for $r>r_0$. All that is required is to join the two solutions at $r=r_0$ such that they are continuous, and their derivative is discontinuous with a jump exactly matching the amplitude of the delta function. This may be written succinctly as

(3.1)\begin{equation} \tilde{G}=\frac{\omega-U(r_0)\,k}{2{\rm \pi}\mathrm{i} r_0}\,\frac{\psi_1(\check{r})\, \psi_2(\hat{r})}{\mathcal{W}(\psi_1,\psi_2;r_0)}, \end{equation}

where

(3.2a,b)\begin{equation} \hat{r}=\max(r,r_0), \quad \check{r}=\min(r,r_0), \end{equation}

and once again $\mathcal {W}(\psi _1,\psi _2;r)$ is the Wronskian of $\psi _1$ and $\psi _2$. Using (2.29), this may be rewritten as

(3.3)\begin{equation} \tilde{G}=\frac{\omega-U(r^{*})\,k}{2{\rm \pi}\mathrm{i} r^{*}\,W(r^{*})}\, \frac{\psi_1(\check{r})\,\psi_2(\hat{r})}{C_1\hat{D}_2-\hat{C}_2D_1}, \quad \text{where } r^{*}=\max(1-h,r_0). \end{equation}

3.2. Analytic continuation behind the critical layer branch cut

The solution for $\tilde {G}$ in (3.3) contains a branch cut along the critical layer $k\in [{{\omega }/{M}},\infty )$. We now introduce the following additional notation. When evaluating a function $f(k)$ on the branch cut, for $k\in [{{\omega }/{M}},\infty )$, we denote

(3.4a–c)\begin{equation} f^{+}(k) = \lim_{\varepsilon\to 0}f(k+\mathrm{i}\varepsilon), \quad f^{-}(k) = \lim_{\varepsilon\to 0}f(k-\mathrm{i}\varepsilon),\quad {\rm \Delta} f(k) = f^{+}(k) - f^{-}(k). \end{equation}

Note that the definition of ${\rm \Delta} f$ agrees with the use of ${\rm \Delta}$ in (2.17), (2.24) and (2.25). By using these equations, we find that

(3.5)\begin{align} {\rm \Delta} \tilde{G}&={-}\frac{\omega-U(r^{*})\,k}{2\mathrm{i}{\rm \pi} r^{*}\,W(r^{*})}\, \frac{1}{C^{-}_1\hat{D}_2-\hat{C}_2D_1+2\mathrm{i}{\rm \pi} AD_1\hat{D}_2} \nonumber\\ &\quad\times \left[\frac{2\mathrm{i}{\rm \pi} A D_1\hat{D}_2\,\psi^{-}_1(\check{r})\, \psi^{-}_2(\hat{r})}{C_1^{-}\hat{D}_2-\hat{C}_2D_1}-\psi^{-}_1(\check{r})\,{\rm \Delta} \psi_2(\hat{r})-{\rm \Delta} \psi_1(\check{r})\,\psi^{-}_2(\hat{r})-{\rm \Delta} \psi_1(\check{r})\,{\rm \Delta} \psi_2(\hat{r})\right]. \end{align}

A typical branch cut, such as the branch cut in $\sqrt {z-z_0}$, may be taken in any direction from the branch point $z_0$. The critical layer branch cut in the complex $k$-plane is different, in that the choice of branch cut was forced upon us by the requirement that the solution be continuous in $r$ for $r\in [1-h,1]$. Nonetheless, noting from (3.5) that ${\rm \Delta} \tilde {G}$ is a well-defined function for general complex $k$, we may use (3.5) to analytically continue $\tilde {G}$ behind the critical layer branch cut. For real $\omega$, we therefore define the analytic continuation of $\tilde {G}$ behind the branch cut into the lower half $k$-plane as

(3.6)\begin{equation} \tilde{G}^{+}(k) = \begin{cases} \tilde{G}(k), & \displaystyle {\rm Im} (k) > 0 \text{ or } {\rm Re}(k) < \dfrac{\omega}{M},\\ \tilde{G}(k)+{\rm \Delta} G(k), & \displaystyle {\rm Im} (k) < 0 \text{ and } {\rm Re}(k) > \dfrac{\omega}{M}. \end{cases} \end{equation}

Similarly, we may rewrite (3.5) as

(3.7)\begin{align} {\rm \Delta} \tilde{G}&={-}\frac{\omega-U(r^{*})\,k}{2\mathrm{i}{\rm \pi} r^{*}\,W(r^{*})}\, \frac{1}{C^{+}_1\hat{D}_2-\hat{C}_2D_1-2\mathrm{i}{\rm \pi} AD_1\hat{D}_2}\nonumber\\ &\quad\times\left[\frac{2\mathrm{i}{\rm \pi} A D_1\hat{D}_2\,\psi^{+}_1(\check{r})\, \psi^{+}_2(\hat{r})}{C_1^{+}\hat{D}_2-\hat{C}_2D_1}-\psi^{+}_1(\check{r})\,{\rm \Delta} \psi_2(\hat{r})-{\rm \Delta} \psi_1(\check{r})\,\psi^{+}_2(\hat{r})+{\rm \Delta} \psi_1(\check{r})\,{\rm \Delta} \psi_2(\hat{r})\right], \end{align}

which allows the analytic continuation of $\tilde {G}$ into the upper half $k$-plane,

(3.8)\begin{equation} \tilde{G}^{-}(k) = \begin{cases} \tilde{G}(k), & \displaystyle {\rm Im} (k) < 0 \text{ or } {\rm Re}(k) < \dfrac{\omega}{M},\\ \tilde{G}(k)-{\rm \Delta} G(k), & \displaystyle {\rm Im} (k) > 0 \text{ and } {\rm Re}(k) > \dfrac{\omega}{M}. \end{cases} \end{equation}

The utility of these analytic continuations is not readily apparent. However, their use allows for poles of $\tilde {G}$, corresponding to modal solutions of the homogeneous Pridmore-Brown equation, to be tracked behind the branch cut, and in particular, a possible hydrodynamic instability mode will later be found to be hidden behind the critical layer branch cut in certain cases. Their use also allows the deformation of integral contours behind the critical layer branch cut, as will be needed for the steepest descent contours needed for the large-$x$ asymptotic evaluation of the inverse Fourier transform.

In what follows, $k^{+}$ and $k^{-}$ denote modal poles (see § 2.5) of only $\tilde {G}^{+}$ or $\tilde {G}^{-}$ respectively.

3.3. Inverting the Fourier transform

Having formulated $\tilde {G}$ as the solution of the inhomogeneous Pridmore-Brown equation (2.5a), to recover the actual pressure perturbation $\hat {p}(x,r,\theta )$, we are required to invert the Fourier transform and sum the Fourier series. For a fixed azimuthal mode $m$, we invert the Fourier transform using the formula

(3.9)\begin{equation} G(x,r;r_0,m)=\frac{1}{2{\rm \pi}}\int_{\mathcal{C}} \tilde{G}(r,r_0,k,m)\, \mathrm{e}^{-\mathrm{i} k x}\,{\rm d}k. \end{equation}

Note, however, that the critical layer branch cut is located along the real $k$-axis, $k\in [{{\omega }/{M}},\infty )$. We are therefore required to be careful in choosing a suitable inversion contour $\mathcal {C}$.

3.3.1. Choosing an inversion contour

In order to choose the correct Fourier inversion contour $\mathcal {C}$, we appeal to the Briggs–Bers criterion (Briggs Reference Briggs1964; Bers Reference Bers1983). The Briggs–Bers criterion, summarized below, invokes the notion of causality: that the cause of the disturbance (the delta function forcing) should occur before the effect (the disturbance $\hat {p}$), which is otherwise lost when considering a time-harmonic forcing, as we do here. A more in-depth description is available in many places in the literature (e.g. Brambley Reference Brambley2009, Appendix A).

In order to make use of the Briggs–Bers criterion, the rate of exponential growth of the solution must be bounded; that is, there must exist $\varOmega,K>0$ such that if $\textrm {Im} (\omega )<-\varOmega$, then $\tilde {G}$ is analytic for $|\textrm {Im} (k)|< K$. For a given $\omega$ with $\textrm {Im} (\omega )<-\varOmega$, we take the $k$-inversion contour $\mathcal {C}$ in (3.9) along the real $k$-axis, and map the locations of any singularities (e.g. poles, branch points, etc.). In order to find a correct integration contour for the real values of $\omega$ that are of interest, the imaginary part of $\omega$ is increased smoothly to $0$, and the locations of any singularities tracked throughout this process. During this process, the $k$-inversion contour $\mathcal {C}$ must be deformed smoothly in order to maintain analyticity; that is, no singularities must cross the $k$-inversion contour. Assuming that this process may be completed and $\textrm {Im} (\omega )$ increased to zero, the resulting $k$-inversion contour $\mathcal {C}$ is the correct causal contour. Since for $x<0$ the $\exp \{-\mathrm {i} kx\}$ term is exponentially small as $|k| \to \infty$ in the upper half $k$-plane, for $x<0$ we may close the contour with a large semicircular arc at infinity in the upper half $k$-plane, denoted $\mathcal {C}_>$. The resulting contours, for real $\omega$, are illustrated in figure 3 for a typical unstable case. The majority of singularities of $\tilde {G}$ are poles that do not cross the real $k$-axis as $\textrm {Im} (\omega )$ is varied, and hence correspond to exponentially decaying disturbances away from the point mass source at $x=0$. The exception to these poles is the pole labelled $k^{+}$, which for this illustration originates in the lower half $k$-plane for $\textrm {Im} (\omega )$ sufficiently negative, and therefore belongs below the $k$-inversion contour. This implies that this pole is seen downstream of the point mass source, for $x>0$, despite having $\textrm {Im} (k)>0$, and therefore corresponds to an exponentially growing instability. For a typical stable case, the situation is the same as shown in figure 3 but with the $k^{+}$ pole not present. Irrespective of the stability, the critical layer, as described earlier, exists when $k/\omega = 1/U(r_c) \in [1/M, \infty ]$ for some critical radius $r_c$, and so is found in the lower half $k$-plane for $\textrm {Im} (\omega )<0$. Thus, as shown in figure 3, for $x>0$ in order to close $\mathcal {C}$ in the lower half $k$-plane, we must pass around the critical layer branch cut, denoted by the contour $\mathcal {C}_b$, before closing in the lower half $k$-plane with a semicircular arc denoted $\mathcal {C}_>$. The contribution from integrating around the critical layer branch cut $\mathcal {C}_b$ leads to the non-modal contribution of the critical layer, and is discussed in detail in § 3.3.3.

Figure 3.

Illustration of typical pole locations, branch cuts and inversion contours taken when an unstable $k^{+}$ pole is present for real $\omega$. The inversion contour for $\tilde {G}$ is labelled $\mathcal {C}$. (a) For $x<0$, the contour is closed in the upper half-plane along the $\mathcal {C}_<$ contour. (b) For $x>0$, the contour is closed in the lower half-plane along the $\mathcal {C}_>$ contour, and around the critical layer branch cut along the $\mathcal {C}_b$ contour. Contributing modal poles are indicated in blue.

3.3.2. Contribution from the poles of $\tilde {G}$

We may now write the integral around the closed contour as a sum of residues of poles:

(3.10a)\begin{gather} \frac{1}{2{\rm \pi}}\int_{{\mathcal{C}\cup\mathcal{C}_<}} \tilde{G}(r,r_0,k,m)\,\mathrm{e}^{-\mathrm{i} k x}\, {\rm d}k = G(x,r;r_0,m) = \sum_{{j:\ {\rm Im} (k_j)>{\rm Im} (\mathcal{C})}} R(k_j) \quad \text{for } x<0, \end{gather}

(3.10b)\begin{gather}\frac{1}{2{\rm \pi}}\int_{{\mathcal{C}\cup\mathcal{C}_b\cup\mathcal{C}_>}} \tilde{G}(r,r_0,k,m)\,\mathrm{e}^{-\mathrm{i} k x}\, {\rm d}k = G(x,r;r_0,m) - I(x) = \sum_{{j:\ {\rm Im} (k_j) < {\rm Im} (\mathcal{C})}} R(k_j)\quad \text{for } x>0, \end{gather}

where $I(x)$ is the contribution from integrating around the critical layer branch cut contour $\mathcal {C}_b$ discussed in the next subsubsection, $R(k_j)$ is the residue from a pole at $k_j$ discussed below, and the notation $\textrm {Im} (k_j) > \textrm {Im} (\mathcal {C})$ is used to denote poles $k_j$ lying above the inversion contour $\mathcal {C}$.

The poles of $\tilde {G}$ correspond to zeros of the denominator of $\tilde {G}$, as given in (3.3). They can occur in two ways: as modal or non-modal poles. We consider the modal poles first. The modal poles occur as zeros of the term $C_1\hat {D}_2-\hat {C}_2 D_1=0$. As discussed in § 2.5, this occurs when both $\psi _1$ and $\psi _2$ satisfy the boundary conditions at $r=0$ and $r=1$. These modal poles can be classified further into acoustic modes and surface modes: acoustic modes are those for which $\alpha$ in (2.11) has a small imaginary part, and correspond to functions that are oscillatory in $r$; surface modes are those for which $\alpha$ has a significant imaginary part, and correspond to functions that decay exponentially away from the duct walls at $r=1$. For different parameters, we may find a variety of surface modes, and two with which we will be particularly interested here will be denoted $k^{-}$ and $k^{+}$. For further details of surface modes, the reader is referred to the existing literature (e.g. Rienstra Reference Rienstra2003; Brambley Reference Brambley2013).

Since the modal poles occur as zeros of $C_1\hat {D}_2-\hat {C}_2 D_1=0$, which we will assume are simple zeros, the contributions from the residues of these poles are given by

(3.11)\begin{equation} R(k)={-}\mathrm{sgn}(x)\,\frac{\omega-U(r^{*})\,k}{2{\rm \pi} r^{*}\,W(r^{*})}\, \frac{\psi_1(\check{r})\,\psi_2(\hat{r})}{\displaystyle \frac{\partial {} }{\partial {k} } (C_1\hat{D}_2-\hat{C}_2D_1)}\,\mathrm{e}^{-\mathrm{i} k x}. \end{equation}

The second type of poles are the non-modal poles, which occur when $W(r^{*})=0$. These occur when we lose independence between $\tilde {p}_1$ and $\tilde {p}_2$ at $r^{*}$. Note from the formula (2.20) for $W(r)$ that $W(r^{*}) = 0$ implies that $r^{*} = r_c^{+}$. Since $1-h\leq r^{*}\leq 1$, this can occur only when $k$ is located on the critical layer branch cut. In what follows, we will refer to this non-modal pole as $k_0$. Note that $k_0$ is a function of the radial location of the point source $r_0$ (through $r^{*}$), which is unlike the modal poles for which $k_j$ is independent of the value of $r_0$; this is one reason why this $k_0$ pole is referred to as a non-modal pole. However, since our closed contour goes around the critical layer branch cut (along contour $\mathcal {C}_b$), this pole is always excluded from the sum of residues in (3.10), and occurs only within the calculation of $I(x)$, which we consider next.

3.3.3. Contribution from the critical layer branch cut

The contribution from the critical layer branch cut, including any non-modal pole $k_0$ along the branch cut, is contained solely within the integral along the critical layer branch cut denoted $\mathcal {C}_b$ in figure 3:

(3.12)\begin{equation} I(x) = \frac{-1}{2{\rm \pi}}\int_{\mathcal{C}_b} \tilde{G} \,\mathrm{e}^{-\mathrm{i} k x}\,\mathrm{d} k. \end{equation}

However, as it stands, this integral for $I(x)$ is oscillatory, owing to the $\mathrm {e}^{-\mathrm {i} kx}$ factor in the integrand, so it is difficult to accurately compute numerically. This is especially true for large values of $x$. Instead, it is helpful to deform the integral onto the steepest descent contour, for which $\mathrm {e}^{-\mathrm {i} kx}$ is exponentially decaying along the contour. This contour deformation has three benefits: first, it allows accurate numerical calculation of the integral; second, it allows the derivation of large-$x$ asymptotics using the method of steepest descent; and third, it brings insight into the various contributions that make up $I(x)$. In deforming the integration contour, however, we must analytically continue $\tilde {G}$ behind the branch cut (as described in § 3.2), and carefully deform around any poles and branch points. This is illustrated schematically in figure 4. Note that poles and branch points of $\tilde {G}$ may exist behind the critical layer branch cut, and we must therefore use analytic continuations of $\tilde {G}$; the reader is reminded that $\tilde {G}^{+}$ is the analytic continuation of $\tilde {G}$ down behind the branch cut from above, while $\tilde {G}^{-}$ is the analytic continuation of $\tilde {G}$ up behind the branch cut from below. Here, we use the notation that a pole of $\tilde {G}^{+}$ with $\textrm {Re}(k)>{{\omega }/{M}}$ is denoted $k^{+}$, and a pole of $\tilde {G}^{-}$ with $\textrm {Re}(k)>{{\omega }/{M}}$ is denoted $k^{-}$. Thus a $k^{+}$ pole with $\textrm {Im} (k^{+})<0$, or a $k^{-}$ pole with $\textrm {Im} (k^{-})>0$, are considered as being hidden behind the critical layer branch cut. In the schematic in figure 4, one $k^{-}$ and one $k^{+}$ pole are present, both with $\textrm {Im} (k)<0$, although this is not always the case; moreover, if present, the $k^{+}$ pole may interact with the integral contours around $k_<$ and $k_>$, in addition to interacting with the integral contour around ${{\omega }/{M}}$ depicted in figure 4, depending on the location of the $k^{+}$ pole.

Figure 4.

(a) Illustration of the integration contour required for the computation of the contribution from the critical layer branch cut, understood by integrating above and below the branch cut. Possible poles of $\tilde {G}^{-}$ and $\tilde {G}^{+}$ are denoted $k^{-}$ and $k^{+}$, respectively. (b) The integration contour after being transformed onto the steepest descent contour. Red lines behave as if evaluated below ${{\omega }/{M}}$ (using $\tilde {G}^{-}$); blue as if having been analytically continued around the ${{\omega }/{M}}$ branch point; green as if having been analytically continued around the ${{\omega }/{M}}$ and $k_<$ branch points; and purple as if having been analytically continued around all branch points, giving $\tilde {G}^{+}$. Note that we have been required to deform contours around the $k^{+}$ and $k^{-}$ poles.

The steepest descent contours are where $\mathrm {e}^{-\mathrm {i} kx}$ is decaying exponentially, i.e. towards $-\mathrm {i}\infty$ in the complex $k$-plane. There is no difficulty deforming the contour at infinity, since $\mathrm {e}^{-\mathrm {i} kx}$ is exponentially small there (provided that $x>0$, which is the only case in which the critical layer branch cut contributes). Along the branch cut, there are up to three branch points singularities, denoted ${{\omega }/{M}}$, $k_<$ and $k_>$ in figure 4, that must be deformed around. These occur because of the presence of the $\log (r-r_c^{+})$ term in $\tilde {p}_2(r)$, and the presence of $\tilde {p}_2(1 - h)$, $\tilde {p}_2(r_0)$ and $\tilde {p}_2(r)$ in the expression for $\tilde {G}$; each of these terms leads to a branch point, respectively at ${{\omega }/{M}}$, at $k_0$ corresponding to $r_c^{+}(k_0) = r_0$, and at $k_r$ corresponding to $r_c^{+}(k_r) = r$.

Moreover, $\tilde {G}$ possesses a pole at $k_0$, which is exactly the non-modal pole referred to above, although there are no poles of $\tilde {G}$ at ${{\omega }/{M}}$ or at $k_r$. Details of these calculations are given in Appendix C. The branch point at $k_r$ is not present when $r \leq 1-h$, and the pole and branch point at $k_0$ are not present when $r_0 \leq 1-h$. For simplicity in what follows, we denote $k_< = \min \{k_0, k_r\}$ and $k_> = \max \{k_0, k_r\}$, as depicted in figure 4.

The total integral around the branch cut can therefore be found by summing these three integrals, subtracting any $k^{-}$ contributions below the branch cut and adding any $k^{+}$ contributions below the branch cut, and adding the pole residue at $k_0$ calculated as if it was located above the branch cut. This results in

(3.13)\begin{equation} I(x)=I_{{{\omega}/{M}}}(x)+I_0(x)+I_r(x)+R_0^{+}(k_0)+ \sum_{{{\rm Im} (k^{+})<0}}R^{+}(k^{+})-\sum_{{{\rm Im} (k^{-})<0}}R^{-}(k^{-}), \end{equation}

where $R^{\pm }$ is the residue given in (3.11) evaluated using $\tilde {G}^{\pm }$. Here, $R_0^{+}(k_0)$ is the residue of the non-modal pole $k_0$ evaluated as if approached from above the branch cut, derived in § C.2 and given in (C20) as

(3.14)\begin{equation} R_0^{+}(k_0) = \frac{2Mk_0^{2}(\omega-Mk_0)\,\mathrm{e}^{-\mathrm{i} k_0x}}{3{\rm \pi} r_0h^{2}\omega(C_1^{+}\hat{D}_2-\hat{C}_2D_1)} \times\begin{cases} \hat{D}_2\,\psi_1(r), & r < r_0,\\ D_1\,\psi_2(r), & r > r_0. \end{cases} \end{equation}

The steepest descent integrals are defined as

(3.15)\begin{equation} I_q(x)=\frac{1}{2{\rm \pi}\mathrm{i}}\int_0^{\infty} {\rm \Delta} \tilde{G}_q(k_q-\mathrm{i}\xi)\,\mathrm{e}^{-\mathrm{i}(k_q-\mathrm{i}\xi)x}\,\mathrm{d}\xi, \end{equation}

and the jumps across each of the steepest descent branch cuts are calculated in Appendix B to be

(3.16a)\begin{align} {\rm \Delta} \tilde{G}_{{{\omega}/{M}}}&=\frac{-(\omega-U(r^{*})\,k) A }{ r^{*}\,W(r^{*})\, (C_1^{-}\hat{D}_2-\hat{C}_2D_1+2\mathrm{i}{\rm \pi} AD_1\hat{D}_2)}\nonumber\\ &\quad\times\begin{cases} \dfrac{\hat{D}_2^{2}\,\psi_1^{-}(\check{r})\,\psi_1^{-}(\hat{r})}{{C_1^{-}\hat{D}_2-\hat{C}_2D_1}}, & \hat{r}<1 - h, \\ \dfrac{D_1\hat{D}_2\,\psi_1^{-}(\check{r})\,\psi_2^{-}(\hat{r})}{{C_1^{-}\hat{D}_2-\hat{C}_2D_1}}, & \check{r}<1 - h<\hat{r}, \\ \dfrac{D_1^{2}\,\psi_2^{-}(\check{r})\,\psi_2^{-}(\hat{r})}{{C_1^{-}\hat{D}_2-\hat{C}_2D_1}}, & 1 - h<\check{r}, \end{cases} \end{align}

(3.16b)\begin{align}{\rm \Delta} \tilde{G}_< &= \frac{-(\omega-U(r^{*})k)}{ r^{*}\,W(r^{*})}\, \frac{ A D_1\,\tilde{p}_1(\check{r})\,\psi_2^{-}(\hat{r}) }{C_1^{-}\hat{D}_2-\hat{C}_2D_1+ 2\mathrm{i}{\rm \pi} A D_1\hat{D}_2}\,H(\check{r}-(1-h)), \end{align}

(3.16c)\begin{align}{\rm \Delta} \tilde{G}_> &= \frac{-(\omega-U(r^{*})\,k)}{ r^{*}\,W(r^{*})} \frac{ A \hat{D}_2\,\psi_1^{-}(\check{r})\,\tilde{p}_1(\hat{r}) }{C_1^{-}\hat{D}_2-\hat{C}_2D_1+2\mathrm{i} {\rm \pi}A D_1\hat{D}_2}\,H(\hat{r}-(1-h)). \end{align}

Note that ${\rm \Delta} \tilde {G}_{{{\omega }/{M}}} + {\rm \Delta} \tilde {G}_< + {\rm \Delta} \tilde {G}_> = {\rm \Delta} \tilde {G} = \tilde {G}^{+} - \tilde {G}^{-}$. While these integrals are now amenable to numerical integration, additional understanding of the contribution from the three steepest descent contours may be gained by considering the large-$x$ limit.

3.4. Far-field decay rates of the critical layer contribution

The critical layer branch cut contribution (3.13) contains integrals $I_q(x)$ given by (3.15) that are amenable to asymptotic analysis in the limit $x \to \infty$, using the method of steepest descent. Having already deformed the integration contours onto the steepest descent contours, so that the integrands have had their $x$-dependent oscillation removed and are now exponentially decaying along the contour, we may directly apply Watson's lemma (Watson Reference Watson1918). If some function $q(k)$ satisfies $f(k_q-\mathrm {i}\xi )\sim B\xi ^{\nu }+O(\xi ^{\nu +1})$ to leading order for small $\xi$ with $\nu >-1$, then Watson's lemma implies that for large $x$,

(3.17)\begin{equation} \frac{1}{2\mathrm{i}{\rm \pi}}\int_0^{\infty} f(k_q-\mathrm{i}\xi)\, \mathrm{e}^{-\mathrm{i}(k_q-\mathrm{i}\xi)x}\,\mathrm{d}\xi\sim\frac{B\,\varGamma(\nu+1)\, \mathrm{e}^{-\mathrm{i} k_q x}}{2\mathrm{i}{\rm \pi} x^{\nu+1}}+O(x^{-(\nu+2)}), \end{equation}

where $\varGamma$ is the Gamma function, and in particular, $\varGamma (\nu +1) = \nu !$ for integer $\nu$. For each of the $I_q(x)$ integrals, this can then be interpreted as an algebraically decaying wave of phase velocity ${\omega }/{k_q}$.

In order to find the decay rates of the steepest descent contours, we are required to understand the behaviour of $\psi _1(r,k_q-\mathrm {i}\xi )$ and $\psi _2(r,k_q-\mathrm {i}\xi )$ for small $\xi$ at ${r \in \{ 1-h,r,r_0,1\}}$. Details of this can be found in Appendix C. The result, given in (C13) and (C14), is that for $k={{\omega }/{M}}-\mathrm {i}\xi$, as $\xi \to 0$ with $\xi >0$, we find that

(3.18)\begin{equation} {\rm \Delta} \tilde{G}_{{{\omega}/{M}}}\sim\begin{cases} \xi^{3/2}, & r_0\leq 1-h, \\ \xi^{5/3}, & r_0>1-h. \end{cases} \end{equation}

By Watson's lemma, this results in a wave convected with the flow speed $M = U(1-h)$ and decaying algebraically like $x^{-{5}/{2}}$ when the source is within the region of uniform flow, and like $x^{-{7}/{2}}$ for a source located in the sheared flow; the pre-factor in each case is different, and is also governed by the above expressions.

In the case $r_0>1-h$, the leading-order contribution to ${\rm \Delta} \tilde {G}_0$ as $k\to k_0$ is derived in § C.2 as

(3.19)\begin{equation} {\rm \Delta} \tilde{G}_0 = \frac{A\omega h^{2}\,U(r_0)}{6r_0Mk_0^{2}(r_0-1+h)}\, \frac{(k-k_0)^{2}}{C_1^{-}\hat{D}_2-\hat{C}_2D_1+2{\rm \pi}\mathrm{i} AD_1\hat{D}_2} \times\begin{cases} \hat{D}_2\,\psi_1^{-}(r), & r_0 > r,\\ D_1\,\psi_2^{-}(r), & r_0 < r. \end{cases} \end{equation}

By Watson's lemma, this results in a wave convected with the flow speed at the point source, $U(r_0)$, and decaying algebraically like $x^{-3}$.

Finally, considering ${\rm \Delta} \tilde {G}_r$ as $k\to k_r$, it is found in § C.3 that

(3.20) \begin{align} {\rm \Delta} \tilde{G}_r&\sim\frac{A(\omega-U(r^{*})\,k_r)\omega^{3} h^{6}}{8r^{*}\,W(r^{*})\,M^{3}k_r^{6}(r-1+h)^{3} }\,\frac{(k-k_r)^{3}}{C_1^{-}\hat{D}_2-\hat{C}_2D_1+2\mathrm{i}{\rm \pi} A D_1\hat{D}_2}\nonumber\\ &\quad\times \begin{cases} \hat{D}_2\,\psi_1^{-}(r_0), & r_0 < r,\\ D_1\,\psi_2^{-}(r_0), & r_0 > r. \end{cases} \end{align}

By Watson's lemma, this results in a wave convected with the flow speed $U(r)$ and decaying algebraically like $x^{-4}$.

It may be noted that the decay rates for $I_0$ and $I_r$ are the same as predicted for a linear boundary layer flow profile by Brambley et al. (Reference Brambley, Darau and Rienstra2012a). We now proceed to compare these results with the previous literature.

3.4.1. Comparisons with previous far-field scalings

Our results for the large-$x$ decay of the various components of the critical layer are compared to those predicted by Swinbanks (Reference Swinbanks1975) for a general flow profile, and those predicted by Brambley et al. (Reference Brambley, Darau and Rienstra2012a) for a constant-then-linear flow profile, in table 1. The $I_0$ integral gives a wave with phase velocity equal to that of the mean flow at the location of the point mass source, $U(r_0)$, provided that the point mass source is in a region of sheared flow, $r_0>1-h$. It can be observed in table 1 that agreement is seen in all three works for $r_0>1-h$. While Swinbanks (Reference Swinbanks1975) did not consider the other cases in detail, this work finds further agreement for the $I_r$ contribution with Brambley et al. (Reference Brambley, Darau and Rienstra2012a). In both the linear and quadratic shear flow cases, when the source is located within the region of sheared flow, the $I_0$ contribution is the slowest decaying term. When the source is located within the uniform flow region, the $I_{{{\omega }/{M}}}$ contribution is the slowest decaying term, although this is matched by the $I_r$ contribution for linear shear. It should noted, however, that when $r_0>1-h$, we have in addition the contribution of the non-modal $k_0$ pole, which does not decay.

Table 1.

Comparison of the different decay rates given for a general flow profile by Swinbanks (Reference Swinbanks1975) and for a linear boundary layer flow profile by Brambley et al. (Reference Brambley, Darau and Rienstra2012a) against those given here for a quadratic boundary layer flow profile.

The difference in the behaviour of the $I_{{{\omega }/{M}}}$ integrals may be understood as having two causes. The first is the difference in behaviour of the constant $A$, given in general as

(3.21)\begin{equation} A={-}\frac{1}{3}\left( {\frac{\omega^{2}}{M^{2}}+\frac{m^{2}}{r_c^{2}}} \right) \left( {\frac{U^{\prime\prime}(r_c)}{U^{\prime}(r_c)}-\frac{1}{r_c}} \right)-\frac{2m^{2}}{3r_c^{3}}. \end{equation}

In the case of linear shear flow, the $U^{\prime \prime }$ term is zero for $k\not ={{\omega }/{M}}$ and the resulting expression is $O(1)$ as $k\to {{\omega }/{M}}$. In the case of a quadratic shear boundary layer, $U^{\prime \prime }$ is non-zero and dominates $A$ as $k\to {{\omega }/{M}}$, providing a factor $(k-{{\omega }/{M}})^{-{1}/{2}}$. The remainder of the differences between the decay rates is explained, for $I_{{{\omega }/{M}}}$, by the fact that $(r_c^{+}-(1-h))\sim (k-{{\omega }/{M}})^{{1}/{2}}$ in the quadratic shear case, whereas for linear shear, $(r_c-(1-h))\sim (k-{{\omega }/{M}})$. For the $I_0$ and $I_r$ contributions, where we do not have $r_c^{+}\to 1-h$, all the other terms are equivalent between the linear and quadratic cases, therefore giving the same eventual asymptotics scalings, although the pre-factors may vary significantly. Further details are given in Appendix D.

4.1. Pole locations

The locations of the poles in the complex $k$-plane for parameter sets A1 and A2 are plotted in figure 5. In addition to the usual acoustic modes (denoted as $\ast$ in figure 5), one $k^{+}$ and one $k^{-}$ pole is found for each parameter set. For parameter set A1, both the $k^{+}$ and $k^{-}$ poles are behind the critical layer branch cut, and so would not be found using conventional numerical methods, although the $k^{+}$ pole does still contribute to the total pressure field through the critical layer branch cut contribution, as described in § 3.3.3. In contrast, for parameter set A2, the $k^{+}$ pole is not behind the branch cut and takes the form of a standard modal pole, in this case a hydrodynamic instability surface wave. The stability of the modal poles is verified from the movement of the poles in the $k$-plane as $\textrm {Im} (\omega )$ is decreased from zero, following the Briggs–Bers criterion (shown in figures 5b,d); note that the critical layer branch cut also moves as a function of $\textrm {Im} (\omega )$. Of particular note is that the $k^{+}$ pole for parameter set A1 emerges from behind the critical layer branch cut as $\textrm {Im} (\omega )$ is reduced from zero, becoming a standard modal pole provided that $\textrm {Im} (\omega )$ is taken sufficiently negative.

Figure 5.

Locations of the poles in the complex $k$-plane for parameter sets A1 (a,b) and A2 (c,d). (a,c) For real $\omega$: acoustic modes with $\textrm {Re}(k)<{{\omega }/{M}}$ ($\ast$); $k^{+}$ poles ($+$); $k^{-}$ poles ($\times$); the critical layer branch cut (solid line); and branch points ${{\omega }/{M}}$ and $k_0$ for $r_0=1-{9h}/{10}$ ($\circ$). (b,d) Trajectories of poles for $-50 < \textrm {Im} (\omega ) < 0$. Poles coloured red (a,c) and solid lines (b,d) denote poles contributing to the modal sum. Poles coloured blue (a,c) and dashed lines (b,d) denote poles hidden behind the branch cut (which varies with $\textrm {Im} (\omega )$) and do not contribute.

As discussed in § 3.3, when the $k^{+}$ pole is located above the branch cut, it is unstable, with a contribution growing exponentially in $x$. When the $k^{+}$ pole is located below the branch cut, we do not see its contribution to the modal sum directly, but instead it contributes as part of the branch cut integral, as seen when deforming onto the steepest descent contour. In this latter case, we would observe a contribution that decays in $x$. In both examples, the $k^{-}$ pole is located above the branch cut and does not contribute towards the Fourier inversion. In the event that this $k^{-}$ pole were located below the branch cut, its contribution would almost exactly cancel the critical layer branch cut contribution, again seen by deforming onto the steepest descent contour; however, the $k^{-}$ pole has not been found below the branch cut for any parameters considered here, unlike in the linear boundary layer profile case, which is investigated further in § 4.3.

Also plotted in figure 5 is the critical layer branch cut for $k > {{\omega }/{M}}$, together with the non-modal $k_0$ pole, which is present only for a point mass source within the boundary layer, $r_0>1-h$. The effects of these non-modal contributions are illustrated in the next subsection.

4.2. Branch cut contributions

Figure 6 illustrates, for three different parameter sets (left to right columns), the differences between the three types of contributions occurring due to the presence of the critical layer branch cut: the sum of the three steepest descent contour integrals (row (i)); the $k_0$ non-modal pole (row (ii)); and the $k^{+}$ modal pole (row (iii)), which is located below the branch cut for all three parameter sets and therefore does not appear in the modal sum. The sum of these contributions is plotted in row (iv) of figure 6. Comparing the non-modal $k_0$ pole (row (ii)) to the sum of the three steepest descent integrals (row (i)), the non-modal $k_0$ pole appears in all three parameter sets to have a small contribution compared to the sum of the three steepest descent integrals for small $x$, although it is comparable and even dominant for larger $x$. For small $x$, the $k^{+}$ pole's contribution (row (iii)) is greater than those of the three steepest descent integrals (row (i)) and the non-modal $k_0$ pole (row (iii)), particularly near the wall at $r=1$. However, since the $k^{+}$ pole decays exponentially in $x$, the non-modal pole will dominate the far-field behaviour of the critical layer branch cut, as can be seen by comparing row (ii) with row (iv).

Figure 6.

A comparison of the terms that contribute to the critical layer, for $r_0=1-{4h}/{5}$. Plotted are the absolute values on a $\log _{10}$ scale. Left to right: parameter sets A1, B and C. Top to bottom: (i) the sum of the three steepest descent contours, $I_{{{\omega }/{M}}}+I_{r}+I_0$; (ii) the non-modal $k_0$ pole; (iii) the contribution of the $k^{+}$ pole located behind the branch cut; and (iv) the total contribution from integrating around the critical layer branch cut, obtained by summing (i)–(iii).

We can look further at how these contributions vary as we adjust the location of the source, shown in figure 7. The contribution of the non-modal $k_0$ pole is seen to be far smaller for the case when $r_0$ is closest to $1-h$ (figure 7a). Note that there is no non-modal $k_0$ pole when $r_0 \leq 1-h$, in which case the dominant contribution will be from the $k^{+}$ pole and the steepest descent contours, which in all cases decay to zero.

Figure 7.

Plots of the real part of the contribution from integrating around the branch cut ($\textrm {Re}(p(x,r))$) for parameter set C, excluding any $k^{+}$ pole located below the branch cut. Solid lines indicate positive values, dashed lines indicate negative values: (a) $r_0=1-{9h}/{10}$; (b) $r_0=1-{3h}/{5}$.

Figure 8 compares the numerically computed steepest descent integrals with their predicted far-field rates of decay given in § 3.4, and good agreement is seen in all cases.

Figure 8.

Plotted for parameter set A1 are $|I_{{{\omega }/{M}}}|$ (top), $|I_r|$ (middle), and $|I_0|$ (bottom). The point source is at $r_0=1-{4h}/{5}$ (a) and $r_0=0.4$ (b). Solid lines correspond to radial locations $r=0.2$, $0.6$, $1-{9h}/{10}$ and $1-{3h}/{5}$. The dashed line is the predicted far-field rate of decay according to § 3.4. Note that for $r=0.2$ and $r=0.6$, $I_r$ is identically zero, since the branch point does not exist. Similarly for $r_0=0.4$ and $I_0$.

4.3. Full Fourier inversion

We now consider the full Fourier inversion, including the contribution from all the modal poles as well as the critical layer branch cut contribution considered above. Figure 9 compares a snapshot in the near field (for small $x$ values) of the wave field generated by only the stable modal poles (figure 9a) with the full solution including the critical layer and any unstable $k^{+}$ pole (figure 9b). When the $k^{+}$ pole is a convective instability, it clearly dominates the solution sufficiently far downstream, as it grows exponentially in $x$. In these near-field plots, the critical layer often appears negligible compared with the modal sum, although in some circumstances it can have a significant effect, as shown by the plots of case C.

Figure 9.

Plotting the real values of the different contributions. (a) Just the contribution for the stable modal poles. (b) The full Fourier inversion, which also includes the $k^{+}$ pole. The parameter sets used from top to bottom are A1, A2, B and C, with $r_0=1-{4h}/{5}$ in each case. In case A2, the $k^{+}$ pole is a convective instability. In cases A1, B and C, the $k^{+}$ pole is located behind the branch cut.

In comparison, figure 10 shows the behaviour outside the near field for the three stable cases from figure 9, plotting the amplitude of oscillations $|p|$ on a logarithmic scale. In all cases, since the modal sum decays exponentially, in the far field, the dominant contribution is from the critical layer, and this is often true from only one or two radii downstream of the point source.

Figure 10.

The absolute value of pressure on a log scale ($10\log _{10}(|p|$)) over a longer range of axial distances downstream of the point source. (a) Just the contribution for the modal poles. (b) Modal poles plus the three steepest descent contours and the $k_0$ non-modal pole. (c) The full Fourier inversion, which also includes the $k^{+}$ pole. The parameter sets used from top to bottom are A1, B and C, with $r_0=1-{4h}/{5}$ in each case. In each case, the $k^{+}$ pole is located behind the branch cut. In the bottom left plot, the contribution from the modal poles is too small to be shown (with $10\log _{10}(|p|) < -78$).

Figure 11 compares the wave field generated in a quadratic boundary layer with the wave field in a linear boundary layer profile (as studied by Brambley et al. Reference Brambley, Darau and Rienstra2012a). The wave fields are reasonably similar, although when the point mass source is within the boundary layer, significant differences are seen downstream. This is related to whether the $k^{+}$ pole is located above or below the branch cut. In the quadratic case, the $k^{+}$ pole always contributes, whether it is behind the branch cut or not, while the $k^{-}$ is always found above the branch cut and so is not seen to contribute at all. With the linear boundary layer, instead we find a $k^{-}$ pole that can be located either above or below the branch cut, while the $k^{+}$ pole is instead located above in all cases. The result of this is that the linear boundary layer profile is always found to be convectively unstable, while the quadratic boundary layer profile is found to be unstable only if the boundary layer is sufficiently thin. Even when both flow profiles give rise to a convective instability, we can see in figure 11 that the growth rates of the instabilities can be significantly different.

Figure 11.

Plotting the real values of the full solution for a quadratic boundary layer flow profile (a) and a linear boundary layer flow profile (b) (from Brambley et al. Reference Brambley, Darau and Rienstra2012a). From left to right, parameters are: set A1 with $r_0=0.4$; set A1 with $r_0=1-{4h}/{5}=0.96$; set A2 with $r_0=0.4$; and set A2 with $r_0=1-{4h}/{5}=0.9992$.

The change in nature of the $k^{+}$ pole in the quadratic case from convective instability to stable is clearly of significant importance. We therefore finally consider the variation of the solution as various parameters of interest are varied.

4.4. Variation of results with changing parameters

The variation of the acoustic modal sum is relatively well understood, so in this subsection we concentrate on the variation of the $k^{+}$ and $k^{-}$ modal poles as various parameters are varied. This includes whether $\textrm {Im} (k^{+})>0$, corresponding to a convective instability, or $\textrm {Im} (k^{+})<0$, corresponding to a stable modal pole hidden behind the branch cut that nonetheless contributes to the modal sum through the branch cut contribution. We also consider whether $\textrm {Im} (k^{-})>0$, meaning the pole is not included, or $\textrm {Im} (k^{-})<0$, in which case the pole is included as part of the contribution of the critical layer branch cut.

Figure 12 illustrates how the $k^{+}$ and $k^{-}$ modal poles vary with boundary layer thickness, frequency, impedance and Mach number. In particular, taking wider boundary layers and lower mean flow velocities appears to stabilize the $k^{+}$ pole, moving it to below the branch cut. In contrast, thinner boundary layers and higher mean flow velocities lead to convective instability. The value of the impedance is also seen to alter the stability of the solution, with, in this case, a range of values of $\textrm {Im} (Z)$ being unstable, and nearly hard-walled values of $|\textrm {Im} (Z)|\to \infty$ being stable, as is seen from the $k^{+}$ poles movement in figure 12(b). The variation of stability as the frequency is varied remains unclear, although it appears likely from figure 12(c) that for certain parameters, there would be a finite range of frequencies for which the $k^{+}$ pole would be unstable, while for frequencies either higher or lower than this range, the $k^{+}$ pole would be stable. Note also from figure 12 that in all cases, the $k^{-}$ pole is located above the branch cut and so does not contribute to either the modal sum or the critical layer branch cut.

Figure 12.

Motion of the modal poles for parameter set C as one parameter is varied (arrows show the motion as the parameter is increased): (a) varying $h$ in $(0.001,0.5)$; (b) varying $\textrm {Im} (Z)$ in $(-\infty,\infty )$, with a dot showing hard-walled values; (c) varying $\textrm {Re}(\omega )$ in $(1,50)$; and (d) varying $M$ in $(0.06,0.9)$. Modal positions for parameter set C are denoted $+$ ($k^{+}$) and $\times$ ($k^{-}$). Dashed lines denote a pole hidden behind the branch cut. Note that (c) and (d) use a rescaled $k$-plane in order for the branch cut to remain fixed as $\omega$ or $M$ is varied.