# The critical layer in quadratic flow boundary layers over acoustic linings

Matthew J. King1  Edward J. Brambley1,2
<EMAIL_ADDRESS>Renan Liupekevicius2,3  Miren Radia4  Paul
Lafourcade5,6  Tauqeer H. Shah7,8  1Mathematics Institute, University of
Warwick, Coventry CV4 7AL, UK 2WMG, University of Warwick, Coventry CV4 7AL,
UK 3Mechanical Engineering, Eindhoven University of Technology, 5600 MB
Eindhoven, NL 4DAMTP, University of Cambridge, Cambridge CB3 0WA, UK 5CEA DAM
Île-de-France, 91297 Arpajon, France 6LMCE, CEA Paris-Saclay, 91680 Bruyères-
le-Châtel, France 7Department of Physics, Government College University
Faisalabad, Faisalabad, Pakistan 8Faculty of Technology, Linnaeus University,
351 95 Växjö, Sweden

###### Abstract

A straight cylindrical duct is considered containing an axial mean flow that
is uniform everywhere except within a boundary layer near the wall, which need
not be thin. Within this boundary layer the mean flow varies parabolically.
The linearized Euler equations are Fourier transformed to give the Pridmore-
Brown equation, for which the Greens function is constructed using Frobenius
series. Inverting the spatial Fourier transform, the critical layer
contribution is given as the non-modal contribution from integrating around
the continuous spectrum branch cut. This contribution is found to be the
dominant downstream contribution to the pressure perturbation in certain
cases, particularly for thicker boundary layers. Moreover, the continuous
spectrum branch cut is found to be involved in stabilizing what are otherwise
convectively unstable modes by hiding them behind the branch cut, particularly
for slower flows. Overall, the contribution from the critical layer is found
to give a neutrally stable non-modal wave with a phase velocity equal to the
mean flow velocity at the source when the source is located within the
sheared-flow region, and to decay algebraically along the duct as
$O(x^{-\frac{5}{2}})$ for a source located with the uniform flow region. The
Frobenius expansion, in addition to being numerically accurate close to the
critical layer where other numerical methods loose accuracy, is also able to
locate modal poles hidden behind the branch cut, which other methods are
unable to find. Matlab code is provided to compute the Greens function.

##  1 Introduction

The propagation of sound through an otherwise steady mean flow has many
important applications. One such application is predicting and optimizing
aircraft engines noise. With aircraft noise being subjected to ever increasing
restrictions, being able to successfully model this noise becomes increasingly
important. In particular, aircraft engine noise at takeoff depends critically
on the sound absorbing performance of acoustic liners. Unfortunately, acoustic
liner performance in the presence of a steady mean flow is poorly predicted by
existing theory, as demonstrated by comparisons to laboratory experiments
(Renou & Aurégan, 2011; Spillere et al., 2020). The theory is equally
applicable to any situation with small perturbations to an otherwise steady
mean flow along a non-rigid boundary: for example, the stability analysis of
flow over a deformable surface.

The behaviour of sound in an otherwise steady mean flow is usually modelled
using the linearized Euler Equations. Non-rigid boundaries, such as the
acoustic liners used in aircraft engines, are usually modelled using an
impedance boundary condition, where a disturbance with oscillating pressure
$\mathrm{Re}(p\exp\\{\mathrm{i}\omega t\\})$ leads to an oscillating normal
boundary velocity $\mathrm{Re}(v\exp\\{\mathrm{i}\omega t\\})$ given by
$p=Z(\omega)v$. Such impedance boundary conditions are well understood for a
mean flow that satisfies no-slip at the boundary. Often, however, we use a
simplified model where the mean flow does not satisfy no-slip at the boundary:
for example, uniform axial flow in a duct. For slipping mean flows, it is
known that the impedance boundary condition must be modified. A common
modified boundary condition is the Myers, or Ingard–Myers, boundary condition
(Ingard, 1959; Myers, 1980). This boundary condition is known to be the
correct limiting behaviour for an inviscid mean flow boundary layer in the
limit that the boundary layer thickness tends to zero (Eversman & Beckemeyer,
1972; Tester, 1973). However, this boundary condition, when applied in the
time domain, is ill-posed (Brambley, 2009). Several alternative boundary
conditions have been suggested (Brambley, 2011b; Schulz et al., 2017; Khamis &
Brambley, 2017; Aurégan, 2018), which each attempt to include more relevant
physics, including the effect of the mean flow boundary layer and the effect
of viscosity. However, these boundary conditions come with their own
complications, including the need to fit further free parameters, and as yet
none have been made to agree with laboratory experiments (Spillere et al.,
2020).

In light of this difficulty with boundary conditions in slipping mean flow,
one may instead only consider mean flows $U(r)$ that satisfy no-slip at the
boundary (e.g. Weng et al., 2017). Doing so, however, involves solving for the
sound in a strongly varying mean flow, which is especially taxing when the
boundary layers are particularly thin. Numerically resolving the sound in thin
boundary layers requires a fine resolution, which then also requires a small
timestep owing to the CFL condition. Progress may be made analytically by
considering the simplified situation of a straight rectilinear or cylindrical
duct containing axial mean flow (as depicted later in figure 1). By
linearizing the Euler Equations about this steady mean flow and assuming
$\exp\\{\mathrm{i}\omega t-\mathrm{i}kx\\}$ dependence, one eventually arrives
at the Pridmore-Brown equation (5), a second-order linear ODE for the pressure
perturbation within the duct due to Pridmore-Brown (1958). The Pridmore-Brown
equation has been the subjected of much analysis (e.g. Mungur & Gladwell,
1969; Ko, 1972; Swinbanks, 1975; Nagel & Brand, 1982; Brambley et al., 2012a;
Rienstra, 2020), owing to its complexity. One complexity is that, treating the
frequency $\omega$ as known and solving for the axial wavenumber $k$ as the
eigenvalue, the Pridmore-Brown equation is not Sturm–Liouville and results in
a nonlinear eigenvalue problem for $k$. A second complexity is that the
Pridmore-Brown equation possesses a regular singularity, referred to as a
_critical layer_ or _continuous spectrum_. Despite these difficulties,
eigenfunction expansions using eigenfunctions of the Pridmore-Brown equation
are frequently used, with the eigenfunctions assumed to form a complete basis
(despite the problem being non-self-adjoint) and the effect of the critical
layer ignored (e.g. Brooks & McAlpine, 2007; Olivieri et al., 2010; Oppeneer
et al., 2016; Rienstra, 2021).

The lack of completeness of the modal solutions of the Pridmore-Brown equation
motivates the investigation of the Green’s function solution. The Green’s
function is the solution of the governing equations subject to a point
forcing; for example, a point mass source leads to the right-hand-side of
equation (5). The Green’s function may be used to construct the solution of
the governing equations subject to any arbitrary forcing; hence, the Green’s
function is capable of being used to express any solution to the governing
equations, in contrast to a modal eigenvalue expansion which can only express
an arbitrary solution if the modal basis is complete. The Green’s function is
also worth considering on its own merits without reference to a particular
forcing, since if the governing equations are capable of exhibiting a
particular feature (such as instability, focusing, perfect reflection, etc),
then the Green’s function must also exhibit that feature. The Green’s function
is also used in various approximation techniques (e.g. Brambley et al., 2012b;
Posson & Peake, 2013; Mathews & Peake, 2018b). For this reason, the Green’s
functions has been constructed for a variety of acoustical situations (e.g.
Rienstra & Tester, 2008; Brambley et al., 2012a; Mathews & Peake, 2017,
2018a). In particular, the Green’s function solution to the Pridmore-Brown
equation naturally includes the critical layer.

The critical layer, or continuous spectrum, is a singularity of the linearized
Euler equations occurring when the phase velocity of the perturbation,
$\omega/k$, is equal to the local fluid velocity of the steady flow,
$U(r_{c})$, for some critical radius $r_{c}$. Because the phase speed is equal
to the flow speed, the effect of the critical layer may be thought of as being
convected with the mean flow, and therefore as hydrodynamic in nature (Case,
1960; Rienstra et al., 2013). For swirling flows, the critical layer is known
to lead to algebraically growing instabilities (Golubev & Atassi, 1996; Tam &
Auriault, 1998; Heaton & Peake, 2006). For the Pridmore-Brown equation, the
critical layer is currently thought to lead to algebraically decaying
disturbances, although publications differ on the exact nature of the decay.
For example, Swinbanks (1975) predicted a disturbance of constant amplitude
plus a disturbance with $O(x^{-3})$ decay for a point source, and $O(x^{-1})$
decay for a distributed source, although exact formulae for these disturbances
are not given. Swinbanks (1975, p. 62) goes on to argue that the constant
amplitude disturbance would not be present when the disturbance is caused “by
moving the surface of a solid body”. In contrast, Félix & Pagneux (2007)
demonstrated numerically, for a point source in a parabolic mean flow, a decay
rate of $O(x^{-1})$. More recently, Brambley et al. (2012a) gave an explicit
analytic solution for the critical layer far-field response for a mean flow
$U(r)$ that is constant in the centre of the duct, and then varies linearly in
a “boundary layer” region to zero at the duct walls. Locating a point source
at a radius $r_{0}$, they found the pressure perturbation from the critical
layer at a radius $r$ consisted of three distinct components with phase
velocities $U(0)$, $U(r)$ and $U(r_{0})$, each with different decay rates.
However, Brambley et al. (2012a) chose a rather special mean flow profile. In
particular, the critical layer is usually caused by a nonzero second
derivative of the mean flow profile, $U^{\prime\prime}(r)$, but for the
constant-then-linear mean flow $U^{\prime\prime}(r)$ is either identically
zero or has a delta function discontinuity; in the constant-then-linear case,
Brambley et al. (2012a) instead attributed the critical layer to the
cylindrical geometry.

In many cases, the effect of the critical layer is negligible in comparison
with the modal sum of the acoustics modes. However, when all acoustic modes
are cut-off and non-propagating, the effect of the critical layer will be
dominant. Moreover, Brambley (2013, figure 6) showed that a mode representing
a hydrodynamic instability could interact with the critical layer, although
this was not seen for a constant-then-linear mean flow profile.

Since the critical layer is a singularity of the Pridmore-Brown equation,
traditional numerical methods are particularly inaccurate near the critical
layer. This often manifests as a collection of spurious numerical modes being
located along the critical layer. In contrast, previous studies have used a
Frobenius expansion about the singular point $r=r_{c}$ (e.g. Heaton & Peake,
2006; Campos & Kobayashi, 2009; Brambley et al., 2012a). This technique both
gives increasing accuracy as the critical layer is approached, and allows
analytical continuation behind the critical layer branch cut. For example,
Brambley et al. (2012a, figure 10) found a previously unknown mode close to
the critical layer that was unable to be resolved numerically using more
traditional finite differences. One complication of the Frobenius series,
however, is that, much like a power series, it has an associated radius of
convergence. For the constant-then-linear mean flow Frobenius expansion
(Brambley et al., 2012a), this did not prove a problem, as the radius of
convergence covered the region of interest in all cases that were considered.
For general flow profiles this will not be the case, and a solution covering
the entire region of interest will involve multiple Frobenius expansions with
overlapping radii of convergence; this will turn out to be the case here. By
matching two different expansions in a region where both converge, a hybrid
solution may be constructed that spans the whole region of interest.

Here, we use the Frobenius expansion method as described by Brambley et al.
(2012a), and apply it to a mean flow that is constant in the centre of the
duct and then varies quadratically within a boundary layer to satisfy non-slip
at the wall. As well as being more realistic than the constant-then-linear
profile considered by Brambley et al. (2012a), this mean flow profile is twice
differentiable, allowing $U^{\prime\prime}(r)$ to enter the analysis, and as
such we expect the results to be more representative of an arbitrary mean flow
profile. The Frobenius expansion is derived in section 2, along with a
derivation of the Pridmore-Brown equations by spatially Fourier transforming
the linearized Euler equations. The Frobenius expansion is then used in
section 3 to derive the Green’s function for a point mass source, including
inverting the spatial Fourier transform and investigating the far-field
behaviour. Results are presented in section 4 by numerically evaluating the
Frobenius expansions and the Green’s function. These results are compared
against previous results, particularly against the predictions by Swinbanks
(1975) and the constant-then-linear results by Brambley et al. (2012a).
Finally, the implications of this work are discussed, and areas for further
research highlighted, in section 5.

##  2 Problem Formulation and Homogeneous Solutions

###  2.1 Constructing the Pridmore-Brown Equation

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

$\displaystyle\frac{\partial\rho}{\partial
t}+\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u})$ $\displaystyle=q,$
$\displaystyle\rho\frac{D\boldsymbol{u}}{Dt}$
$\displaystyle=-\boldsymbol{\nabla}p,$ $\displaystyle\frac{Dp}{Dt}$
$\displaystyle=\frac{\gamma p}{\rho}\frac{D\rho}{Dt}.$ (1)

Here, 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

$q=\mathrm{Re}\left(\frac{\epsilon}{r_{0}}\delta(x)\delta(\theta)\delta(r-r_{0})\exp\\{\mathrm{i}\omega
t\\}\right),$ (2)

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$,

$\displaystyle\rho=\rho_{0}(r)+\mathrm{Re}\big{(}\epsilon\hat{\rho}e^{\mathrm{i}\omega
t}\big{)}+O(\epsilon^{2}),\qquad\qquad
p=p_{0}+\mathrm{Re}\big{(}\epsilon\hat{p}e^{\mathrm{i}\omega
t}\big{)}+O(\epsilon^{2}),$
$\displaystyle\boldsymbol{u}=U(r)\boldsymbol{e_{x}}+\mathrm{Re}\Big{(}\epsilon\big{(}\hat{u},\hat{v},\hat{w}\big{)}e^{\mathrm{i}\omega
t}\Big{)}+O(\epsilon^{2}),$ (3)

where $p_{0}$ is necessarily a constant in order that the steady state should
satisfy the Euler equations. 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

$\hat{p}(x,r,\theta)=\frac{1}{2\upi}\sum_{m=-\infty}^{\infty}\mathrm{e}^{-\mathrm{i}m\theta}\int_{-\infty}^{\infty}\widetilde{p}(r;k,m,\omega)\mathrm{e}^{-\mathrm{i}kx}\,\mathrm{d}k,$
(4)

and similarly for the density $\hat{\rho}$ and the velocity components
$\hat{u}$, $\hat{v}$ and $\hat{w}$. Substituting these into the Euler
equations (1), 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 ODE in the radial
coordinate $r$ for $\widetilde{p}$,

$\displaystyle\widetilde{p}^{\prime\prime}+\left(\frac{2kU^{\prime}}{\omega-U(r)k}+\frac{1}{r}-\frac{\rho^{\prime}_{0}}{\rho_{0}}\right)\\!\widetilde{p}^{\prime}+\left(\frac{(\omega-U(r)k)^{2}}{c_{0}^{2}}-k^{2}-\frac{m^{2}}{r^{2}}\right)\\!\widetilde{p}=\frac{\omega-U(r_{0})k}{2\mathrm{i}\upi
r_{0}}\delta(r-r_{0}),$ (5a)
$\displaystyle\text{with}\qquad\qquad\tilde{v}=\frac{\mathrm{i}\widetilde{p}^{\prime}}{\rho_{0}(\omega-
Uk)},$ (5b)

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

One boundary condition to (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$

$\displaystyle\widetilde{p}(0)$ $\displaystyle=0\quad\text{for }m\not=0$
$\displaystyle\widetilde{p}^{\prime}(0)$ $\displaystyle=0\quad\text{for }m=0$
(6)

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$,

$\displaystyle\widetilde{p}(a)$ $\displaystyle=Z(\omega)\tilde{v}(a)$
$\displaystyle\iff$ $\displaystyle\widetilde{p}^{\prime}(a)$
$\displaystyle=\frac{\mathrm{i}\omega}{Z}\widetilde{p}(a),$ (7)

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

In what follows, we make the simplifying assumptions of a perfect gas and a
constant density $\rho_{0}(r)$. This is a homentropic assumption, and implies
that $c_{0}(r)$ is also constant. We may then nondimensionalize 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 nondimensional
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:

$U(r)=\begin{cases}M&0\leq r\leq 1-h\\\ M(1-(1-\frac{1-r}{h})^{2})&1-h\leq
r\leq 1\end{cases}.$ (8)

With the nondimensionalization of velocities by $c_{0}$, $M$ here is the duct
centreline Mach number. This situation is depicted schematically in figure 1.

$Z=p/v$$x$$r$$h$$U(r)$$\rho_{0}(r)$$\theta$ Figure 1: A cross sectional view
of a cylindrical duct with lined walls containing sheared axial flow.
$\rho_{0}(r)$ is the mean flow density (here taken constant), and $U(r)$ is
the mean flow velocity, here taken to be uniform outside a boundary layer of
width $h$. $Z$ is the boundary impedance and defines the boundary condition at
the wall of the duct.

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

$\widetilde{p}^{\prime\prime}+\left(\frac{2kU^{\prime}}{\omega-U(r)k}+\frac{1}{r}\right)\widetilde{p}^{\prime}+\left((\omega-U(r)k)^{2}-k^{2}-\frac{m^{2}}{r^{2}}\right)\widetilde{p}=0.$
(9)

###  2.2 Homogeneous Solutions Within the Region of Uniform Flow

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

$\widetilde{p}^{\prime\prime}+\frac{1}{r}\widetilde{p}^{\prime}+\left((\omega-
Mk)^{2}-k^{2}-\frac{m^{2}}{r^{2}}\right)\widetilde{p}=0.$ (10)

This is Bessel’s equations of order $m$ rescaled by $\alpha$, where

$\alpha^{2}=(\omega-Mk)^{2}-k^{2};$ (11)

it will turn out later that the branch chosen for $\alpha$ does not matter,
although for definiteness one may choose $\mathrm{Re}(\alpha)>0$. Bessel’s
equation has two pairs of linearly independent solutions that we shall 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, 1964). 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 to the homogeneous Pridmore-
Brown equation (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

$\displaystyle r_{c}^{\pm}$ $\displaystyle=1-h\pm Q$ $\displaystyle Q$
$\displaystyle=h\sqrt{1-\frac{\omega}{Mk}}.$ (12)

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

Following Brambley et al. (2012a), we propose a Frobenius expansion (Teschl,
2012) about the regular singularity $r_{c}^{+}$,

$\displaystyle\widetilde{p}(r)$
$\displaystyle=\sum_{n=0}^{\infty}a_{n}(r-r_{c}^{+})^{n+\sigma}$
$\displaystyle\text{with}\qquad a_{0}\neq 0.$ (13)

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

$\displaystyle\widetilde{p}_{1}(r)$
$\displaystyle=\sum_{n=0}^{\infty}a_{n}(r-r_{c}^{+})^{n+3},$ (14a)
$\displaystyle\widetilde{p}_{2}(r)$
$\displaystyle=A\widetilde{p}_{1}(r)\log(r-r_{c}^{+})+\sum_{n=0}^{\infty}b_{n}(r-r_{c}^{+})^{n}.$
(14b)

The coefficients $a_{n}$ and $b_{n}$ are derived in appendix A, where, in
particular, it is found that

$\displaystyle a_{0}$ $\displaystyle=b_{0}=1,$ $\displaystyle b_{1}$
$\displaystyle=0,$ $\displaystyle b_{2}=$
$\displaystyle-\frac{1}{2}\\!\left(\\!k^{2}+\left(\frac{m}{r_{c}^{+}}\right)^{\\!\\!2}\right)\\!,$
$\displaystyle b_{3}$ $\displaystyle=0,$ (15)

and that

$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}},$
(16)

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

Due to the log term in $\widetilde{p}_{2}$ in (14), 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]$, and 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.

$0$$\mathrm{Re}(r)$$\mathrm{Im}(r)$$1-h$$1$$r_{c}^{+}$$r_{c}^{-}$Branch
cut(a)$0$$\mathrm{Re}(r)$$\mathrm{Im}(r)$$1-h$$1$$r_{c}^{+}$$r_{c}^{-}$Branch
cut(b) Figure 2: Schemetic of possible locations of the $r_{c}^{+}$ branch cut
in the complex $r$-plane. (a) A possible choice of branch cut when
$\mathrm{Im}(r_{c}^{+})>0$. (b) The other choice of branch cut is needed when
$\mathrm{Im}(r_{c}^{+})<0$.

When $r_{c}^{+}$ is real and $1-h<r_{c}^{+}<1$, no suitable choice of branch
cut exists, and as a result any solution $\widetilde{p}(r)$ with
$\widetilde{p}(r_{c}^{+})\neq 0$ necessarily has a singular third derivative
at $r_{c}^{+}$. This only occurs for particular values of $k$, however, and we
can map the corresponding values of $k$ in the complex $k$ plane to find they
fall exactly on the half line $[\frac{\omega}{M},\infty)$; this range of
excluded values of $k$ we refer to as the critical layer branch cut. As
$r_{c}^{+}$ becomes real, note that the value of
$\widetilde{p}_{2}(r_{c}^{+})$ is different depending on whether we approach
from positive or negative imaginary part. Thinking of $r_{c}^{+}(k)$ as a
function of $k$, this corresponds to approaching the critical layer branch cut
$[\frac{\omega}{M},\infty)$ in $k$ from above or below. This re-enforces the
consideration of the critical layer appearing as a branch cut in the complex
$k$ plane, taken along the real line from $k=\frac{\omega}{M}$. The change in
$\widetilde{p}_{2}$ when crossing the critical layer branch cut from below to
above is described as

$\Delta\widetilde{p}_{2}(r)\;=\;\lim_{\mathchoice{\hbox
to0.0pt{\hss$\displaystyle{\mathrm{Im}(k)\searrow 0+}$\hss}}{\hbox
to0.0pt{\hss$\textstyle{\mathrm{Im}(k)\searrow 0+}$\hss}}{\hbox
to0.0pt{\hss$\scriptstyle{\mathrm{Im}(k)\searrow 0+}$\hss}}{\hbox
to0.0pt{\hss$\scriptscriptstyle{\mathrm{Im}(k)\searrow
0+}$\hss}}}\;\widetilde{p}_{2}(r)\;-\;\lim_{\mathchoice{\hbox
to0.0pt{\hss$\displaystyle{\mathrm{Im}(k)\nearrow 0-}$\hss}}{\hbox
to0.0pt{\hss$\textstyle{\mathrm{Im}(k)\nearrow 0-}$\hss}}{\hbox
to0.0pt{\hss$\scriptstyle{\mathrm{Im}(k)\nearrow 0-}$\hss}}{\hbox
to0.0pt{\hss$\scriptscriptstyle{\mathrm{Im}(k)\nearrow
0-}$\hss}}}\;\widetilde{p}_{2}(r)\;=\;-2\pi\mathrm{i}A\widetilde{p}_{1}(r)H(r_{c}^{+}-r)$
(17)

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

In order to retrieve this result we need only consider the $\log(r-r_{c}^{+})$
term of $\widetilde{p}_{2}$. Note that $\partial r_{c}^{+}/\partial k>0$ for
real $k$ and real positive $\omega$; hence, if $k$ is nearly real and
$\mathrm{Im}(k)>0$, then $\mathrm{Im}(r_{c}^{+})>0$, and we must take the
branch cut of $\log(r-r_{c}^{+})$ upwards towards $+\mathrm{i}\infty$.
Similarly, for $\mathrm{Im}(k)<0$ then $\mathrm{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 $-\upi$.
When we choose the branch cut into the lower half plane, this correponds to a
complex argument of $\upi$. This difference results in the jump of
$2\upi\mathrm{i}$ given. If we instead consider $r>r_{c}^{+}$, the same
argument is retrieved regardless of which direction we take the branch-cuts,
and 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 (5) across $r\in[0,1]$,
by matching the solutions derived above in sections 2.2 and 2.3. We construct
$\psi_{1}(r)$ to satisfy the boundary condition at $r=0$ (6), and
$\psi_{2}(r)$ satisfies the boundary condition at $r=1$. Therefore, we take

$\psi_{1}(r)=\begin{cases}J_{m}(\alpha r)&0\leq r\leq 1-h\\\
C_{1}\widetilde{p}_{1}(r)+D_{1}\widetilde{p}_{2}(r)&1-h\leq r\leq
1,\end{cases}$ (18)

where the matching coefficients $C_{1}$ and $D_{1}$ ensure $C^{1}$ continuity,
and are given by

$\displaystyle C_{1}$
$\displaystyle=\phantom{-}\frac{J_{m}(\alpha(1-h))\widetilde{p}_{2}^{\prime}(1-h)-\alpha
J_{m}^{\prime}(\alpha(1-h))\widetilde{p}_{2}(1-h)}{W(1-h)},$ (19a)
$\displaystyle D_{1}$
$\displaystyle={-}\frac{J_{m}(\alpha(1-h))\widetilde{p}_{1}^{\prime}(1-h)-\alpha
J_{m}^{\prime}(\alpha(1-h))\widetilde{p}_{1}(1-h)}{W(1-h)},$ (19b)

and $W(r)=\mathcal{W}(\widetilde{p}_{1},\widetilde{p}_{2};r)$ is the Wronskian
of $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$, given in appendix A.4 as

$W(r)=\mathcal{W}(\widetilde{p}_{1},\widetilde{p}_{2};r)=\widetilde{p}_{1}(r)\widetilde{p}_{2}^{\prime}(r)-\widetilde{p}_{2}(r)\widetilde{p}_{1}^{\prime}(r)=-\frac{3}{4}\frac{r_{c}^{+}(r-r_{c}^{+})^{2}(r-r_{c}^{-})^{2}}{rQ^{2}}.$
(20)

Having constructed $\psi_{1}$ to satisfy the boundary condition at $r=0$, we
now proceed to construct $\psi_{2}$ which satisfies the boundary condition (7)
at $r=1$. Writing $\psi_{2}$ in terms of the homogeneous solutions derived
above,

$\psi_{2}(r)=\begin{cases}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2}H_{m}^{(1)}(\alpha
r)+\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2}H_{m}^{(2)}(\alpha
r)&0\leq r\leq 1-h\\\
\widehat{C}_{2}\widetilde{p}_{1}(r)+\widehat{D}_{2}\widetilde{p}_{2}(r)&1-h\leq
r\leq 1,\end{cases}$ (21)

we choose $\widehat{C}_{2}$ and $\widehat{D}_{2}$ to satisfy $\psi_{2}(1)=1$
and $\psi^{\prime}_{2}(1)=-\frac{\mathrm{i}\omega}{Z}$. This forces a non-zero
normalized solution to $\psi_{2}$ which satisfies the boundary condition (7)
at $r=1$, and leads to

$\displaystyle\widehat{C}_{2}$
$\displaystyle=\frac{{\widetilde{p}_{2}^{\prime}(1)}+\frac{\mathrm{i}\omega}{Z}{\widetilde{p}_{2}(1)}}{W(1)},$
$\displaystyle\widehat{D}_{2}$
$\displaystyle=-\frac{{\widetilde{p}_{1}^{\prime}(1)}+\frac{\mathrm{i}\omega}{Z}{\widetilde{p}_{1}(1)}}{W(1)}.$
(22)

The coefficients
$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2}$
and
$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2}$
are chosen such that our solution is $C^{1}$ continuous at $r=1-h$, giving

$\\!\left(\\!\\!\begin{array}[]{c}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2}\\\
\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2}\end{array}\\!\\!\\!\right)\\!=\frac{\mathrm{i}\upi(1\\!-\\!h)}{4}\\!\left(\\!\\!\\!\begin{array}[]{cc}\phantom{-}\alpha
H_{m}^{(2)\prime}\big{(}\alpha(1\\!-\\!h)\big{)}&\\!\\!\\!{-}H_{m}^{(2)}\big{(}\alpha(1\\!-\\!h)\big{)}\\\
{-}\alpha
H_{m}^{(1)\prime}\big{(}\alpha(1\\!-\\!h)\big{)}&\\!\\!\\!\phantom{-}H_{m}^{(1)}\big{(}\alpha(1\\!-\\!h)\big{)}\end{array}\\!\\!\\!\right)\\!\\!\left(\\!\begin{array}[]{cc}\widetilde{p}_{1}(1\\!-\\!h)&\widetilde{p}_{2}(1\\!-\\!h)\\\
\widetilde{p}_{1}^{\prime}(1\\!-\\!h)&\widetilde{p}_{2}^{\prime}(1\\!-\\!h)\end{array}\\!\\!\right)\\!\left(\\!\\!\begin{array}[]{c}\widehat{C}_{2}\\\
\widehat{D}_{2}\end{array}\\!\\!\\!\right)\\!$ (23)

where the factor at the beginning comes from the Wronskian of $H_{m}^{(1)}$
and $H_{m}^{(2)}$ from Abramowitz & Stegun (1964, 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 $\widetilde{p}_{2}(r)$ when $r<r_{c}^{+}$, we
have, provided $r_{c}^{+}<1$,

$\displaystyle\Delta C_{1}=2\mathrm{i}\upi
AD_{1},\qquad\qquad\Delta\widehat{C}_{2}=\Delta
D_{1}=\Delta\widehat{D}_{2}=0,$ (24a)
$\displaystyle\left(\\!\begin{array}[]{c}\Delta\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2}\\\
\Delta\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2}\end{array}\\!\right)=\frac{\upi^{2}(1\\!-\\!h)A\widehat{D}_{2}}{2}\left(\begin{array}[]{cc}\phantom{-}\alpha
H_{m}^{(2)\prime}\big{(}\alpha(1\\!-\\!h)\big{)}&{-}H_{m}^{(2)}\big{(}\alpha(1\\!-\\!h)\big{)}\\\
{-}\alpha
H_{m}^{(1)\prime}\big{(}\alpha(1\\!-\\!h)\big{)}&\phantom{-}H_{m}^{(1)}\big{(}\alpha(1\\!-\\!h)\big{)}\end{array}\right)\left(\begin{array}[]{c}\widetilde{p}_{1}(1\\!-\\!h)\\\
\widetilde{p}_{1}^{\prime}(1\\!-\\!h)\end{array}\right)$ (24f)

resulting in (provided $r_{c}^{+}<1$)

$\displaystyle\Delta\psi_{1}(r)$ $\displaystyle=2\mathrm{i}\upi
AD_{1}\widetilde{p}_{1}H(r-r_{c}^{+}),$ (25a) $\displaystyle\Delta\psi_{2}(r)$
$\displaystyle=\begin{cases}\Delta\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2}H_{m}^{(1)}(\alpha
r)+\Delta\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2}H_{m}^{(2)}(\alpha
r)&0\leq r\leq 1-h\\\ -2\mathrm{i}\upi
A\widetilde{p}_{1}(r)\widehat{D}_{2}H(r_{c}^{+}-r)&1-h\leq r\leq
1.\end{cases}$ (25b)

Note that, if $r_{c}^{+}>1$, then $\Delta\psi_{1}=\Delta\psi_{2}=0$, since the
$\psi_{1}$ and $\psi_{2}$ solutions are uniquely defined 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 (9) are nonzero
solutions $\widetilde{p}(r)$ that satisfy both the boundary conditions at
$r=0$ and at $r=1$ (6,7). In general, satisfying both boundary conditions
would force the solution $\widetilde{p}(r)\equiv 0$, so nonzero solutions
exist only for particular modal eigenvalues $k$ (assuming $\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,
and so their Wronskian $\mathcal{W}(\psi_{1},\psi_{2};r)$ is not identically
zero, where

$\mathcal{W}(\psi_{1},\psi_{2};r)=\psi_{1}(r)\psi_{2}^{\prime}(r)-\psi_{2}(r)\psi_{1}^{\prime}(r).$
(26)

However, if $\widetilde{p}(r)$ is nonzero and satisfies both boundary
conditions at $r=0$ and $r=1$, then
$\widetilde{p}(r)=a\psi_{1}(r)=b\psi_{2}(r)$ for some nonzero coefficients
$a,b$. In other words, a modal solution is one where $\psi_{1}$ and $\psi_{2}$
are linearly dependent, and so $\mathcal{W}(\psi_{1},\psi_{2};r)\equiv 0$.

For $1-h\leq r\leq 1$, substituting $\psi_{1}$ from (18) and $\psi_{2}$ from
(21) into the Wronskian (26) gives

$\mathcal{W}(\psi_{1},\psi_{2};r)=(C_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1})W(r),$
(27)

where $W(r)$ is the Wronskian between $\widetilde{p}_{1}$ and
$\widetilde{p}_{2}$ and is given earlier in (20). Since $\widetilde{p}_{1}$
and $\widetilde{p}_{2}$ were constructed to be linearly independent, we expect
$W(r)$ not to be identically zero, and indeed (20) shows that $W(r)\neq 0$
except at the critical layer $r=r_{c}^{+}$. A modal solution, therefore, is
given by the condition that $C_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}=0$,
which is independent of $r$, and implies that
$C_{1}/D_{1}=\widehat{C}_{2}/\widehat{D}_{2}$ and 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 (26) becomes

$\mathcal{W}(\psi_{1},\psi_{2};r)\;=\;\alpha\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2}\mathcal{W}(J_{m},H_{m}^{(1)};r)+\alpha\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2}\mathcal{W}(J_{m},H_{m}^{(2)};r)\;=\;\alpha(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2}-\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2})\frac{2\mathrm{i}}{\pi
r},$ (28)

where we have made use of the Bessel function identities 9.1.3, 9.1.4 and
9.1.16 from Abramowitz & Stegun (1964). 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$, since $\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

$\mathcal{W}(\psi_{1},\psi_{2};r)=(C_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1})\begin{cases}W(r)&1-h\leq
r\leq 1\\\ W(1-h)\frac{1-h}{r}&0\leq r\leq 1-h,\end{cases}$ (29)

and that a mode corresponds to the dispersion relation
$0=D(k,\omega)=C_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{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 Inhomogeneous Solutions and Inverting the Fourier Transform

###  3.1 Inhomogeneous Solution to the Pridmore-Brown Equation

While previously we have only been solving the homogeneous form (9), our
original problem was to solve the inhomogeneous Pridmore-Brown equation (5a)
subjected to a harmonic point mass source. Due to the right hand side of (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 $\widetilde{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, $\widetilde{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

$\widetilde{G}=\frac{\omega-U(r_{0})k}{2\upi\mathrm{i}r_{0}}\frac{\psi_{1}(\check{r})\psi_{2}(\hat{r})}{\mathcal{W}(\psi_{1},\psi_{2};r_{0})},$
(30)

where

$\displaystyle\hat{r}$ $\displaystyle=\max(r,r_{0}),$ $\displaystyle\check{r}$
$\displaystyle=\min(r,r_{0}),$ (31)

and once again $\mathcal{W}(\psi_{1},\psi_{2};r)$ is the Wronskian of
$\psi_{1}$ and $\psi_{2}$. Using (29), this may be rewritten as

$\displaystyle\widetilde{G}$
$\displaystyle=\frac{\omega-U(r^{*})k}{2\upi\mathrm{i}r^{*}W(r^{*})}\frac{\psi_{1}(\check{r})\psi_{2}(\hat{r})}{C_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}$
$\displaystyle\text{where}\quad r^{*}$ $\displaystyle=\max(1-h,r_{0})$ (32)

###  3.2 Analytic continuation behind the critical layer branch cut

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

$\displaystyle f^{+}(k)$ $\displaystyle=\lim_{\varepsilon\to
0}f(k+\mathrm{i}\varepsilon)$ $\displaystyle f^{-}(k)$
$\displaystyle=\lim_{\varepsilon\to 0}f(k-\mathrm{i}\varepsilon)$
$\displaystyle\Delta f(k)$ $\displaystyle=f^{+}(k)-f^{-}(k)$ (33)

Note that the definition of $\Delta f$ agrees with the use of $\Delta$ in
equations (17,24,25) above. By using these equations, we find that

$\displaystyle\Delta\widetilde{G}=$
$\displaystyle-\frac{\omega-U(r^{*})k}{2\mathrm{i}\pi
r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}$ (34)
$\displaystyle\quad\times\left[\frac{2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}\psi^{-}_{1}(\check{r})\psi^{-}_{2}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}-\psi^{-}_{1}(\check{r})\Delta\psi_{2}(\hat{r})-\Delta\psi_{1}(\check{r})\psi^{-}_{2}(\hat{r})-\Delta\psi_{1}(\check{r})\Delta\psi_{2}(\hat{r})\right].$

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]$. None-the-less, noting from (34) that $\Delta\widetilde{G}$ is
well defined function for general complex $k$, we may use equation (34) to
analytically continue $\widetilde{G}$ behind the critical layer branch cut.
For real $\omega$, we therefore define the analytic continuation of
$\widetilde{G}$ behind the branch cut into the lower-half $k$-plane as

$\widetilde{G}^{+}(k)=\begin{cases}\widetilde{G}(k)&\mathrm{Im}(k)>0\text{ or
}\mathrm{Re}(k)<\frac{\omega}{M},\\\ \widetilde{G}(k)+\Delta
G(k)&\mathrm{Im}(k)<0\text{ and }\mathrm{Re}(k)>\frac{\omega}{M}.\end{cases}$
(35)

Similarly, we may rewrite (34) as

$\displaystyle\Delta\widetilde{G}=$
$\displaystyle-\frac{\omega-U(r^{*})k}{2\mathrm{i}\pi
r^{*}W(r^{*})}\frac{1}{C^{+}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}-2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}$ (36)
$\displaystyle\quad\times\left[\frac{2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}\psi^{+}_{1}(\check{r})\psi^{+}_{2}(\hat{r})}{C_{1}^{+}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}-\psi^{+}_{1}(\check{r})\Delta\psi_{2}(\hat{r})-\Delta\psi_{1}(\check{r})\psi^{+}_{2}(\hat{r})+\Delta\psi_{1}(\check{r})\Delta\psi_{2}(\hat{r})\right],$

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

$\widetilde{G}^{-}(k)=\begin{cases}\widetilde{G}(k)&\mathrm{Im}(k)<0\text{ or
}\mathrm{Re}(k)<\frac{\omega}{M},\\\ \widetilde{G}(k)-\Delta
G(k)&\mathrm{Im}(k)>0\text{ and }\mathrm{Re}(k)>\frac{\omega}{M}.\end{cases}$
(37)

The utility of these analytic continuations in not readily apparent. However,
their use allows for poles of $\widetilde{G}$, corresponding to modal
solutions to 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 section 2.5, of
only $\widetilde{G}^{+}$ or $\widetilde{G}^{-}$ respectively.

###  3.3 Inverting the Fourier Transform

Having formulated $\widetilde{G}$ as the solution to the inhomogeneous
Pridmore-Brown equation $\eqref{PridmoreBrown}$, 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

$G(x,r;r_{0},m)=\frac{1}{2\upi}\int_{\mathcal{C}}\widetilde{G}(r,r_{0},k,m)\mathrm{e}^{-\mathrm{i}kx}dk.$
(38)

Note, however, that the critical layer branch cut is located along the
real-$k$ axis $k\in[\frac{\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, 1964; Bers, 1983). The Briggs-
Bers criterion, sumarised 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, 2009, appendix A).

In order to make use to the Briggs-Bers criterion, the rate of exponential
growth of the solution must be bounded; that is, there must exist $\Omega,K>0$
such that, if $\mathrm{Im}(\omega)<-\Omega$, then $\widetilde{G}$ is analytic
for $|\mathrm{Im}(k)|<K$. For a given $\omega$ with
$\mathrm{Im}(\omega)<-\Omega$, we take the $k$-inversion contour $\mathcal{C}$
in (38) 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 smoothly increased to $0$, and the locations of any
singularities tracked throughout this process. During this process, the
$k$-inversion contour $\mathcal{C}$ must be smoothly deformed in order to
maintain analyticity; that is, no singularities must cross the $k$-inversion
contour. Assuming this process may be completed and $\mathrm{Im}(\omega)$
increased to zero, then 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 semi-circular arc at infinity in
the upper-half $k$-plane, denoted $\mathcal{C}_{>}$. The resulting contours
are illustrated in 3.

$\frac{\omega}{M}$$\frac{\omega}{M}$$k^{+}$$k^{+}$$\mathcal{C}$$\mathcal{C}$$\mathcal{C_{<}}$$\mathcal{C}_{b}$$\mathcal{C_{>}}$
Figure 3: Illustration of the inversion contours taken when an unstable
$k^{+}$ pole is present. The inversion contour for $\widetilde{G}$ is labelled
$\mathcal{C}$. (Left) For $x<0$, the contour is closed in the upper half plane
along the $\mathcal{C}_{<}$ contour. (Right) 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.

In this illustration, the majority of singularities of $\widetilde{G}$ are
poles which do not cross the real $k$ axis as $\mathrm{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 $\mathrm{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 $\mathrm{Im}(k)>0$, and
therefore corresponds to an exponentially growing instability. 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 $\mathrm{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 semi-circular 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 below in section 3.3.3.

####  3.3.2 Contribution from the poles of $\widetilde{G}$

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

$\displaystyle\frac{1}{2\upi}\int_{\mathchoice{\hbox
to0.0pt{$\displaystyle{\mathcal{C}\cup\mathcal{C}_{<}}$\hss}}{\hbox
to0.0pt{$\textstyle{\mathcal{C}\cup\mathcal{C}_{<}}$\hss}}{\hbox
to0.0pt{$\scriptstyle{\mathcal{C}\cup\mathcal{C}_{<}}$\hss}}{\hbox
to0.0pt{$\scriptscriptstyle{\mathcal{C}\cup\mathcal{C}_{<}}$\hss}}}\;\widetilde{G}(r,r_{0},k,m)\mathrm{e}^{-\mathrm{i}kx}dk$
$\displaystyle\;=\;G(x,r;r_{0},m)\phantom{{}-I(x)}\;=\;\sum_{\mathchoice{\hbox
to0.0pt{\hss$\displaystyle{j:\;\mathrm{Im}(k_{j})>\mathrm{Im}(\mathcal{C})}$\hss}}{\hbox
to0.0pt{\hss$\textstyle{j:\;\mathrm{Im}(k_{j})>\mathrm{Im}(\mathcal{C})}$\hss}}{\hbox
to0.0pt{\hss$\scriptstyle{j:\;\mathrm{Im}(k_{j})>\mathrm{Im}(\mathcal{C})}$\hss}}{\hbox
to0.0pt{\hss$\scriptscriptstyle{j:\;\mathrm{Im}(k_{j})>\mathrm{Im}(\mathcal{C})}$\hss}}}R(k_{j})$
$\displaystyle\text{for}\quad x<0,$ (39a)
$\displaystyle\frac{1}{2\upi}\int_{\mathchoice{\hbox
to0.0pt{$\displaystyle{\mathcal{C}\cup\mathcal{C}_{b}\cup\mathcal{C}_{>}}$\hss}}{\hbox
to0.0pt{$\textstyle{\mathcal{C}\cup\mathcal{C}_{b}\cup\mathcal{C}_{>}}$\hss}}{\hbox
to0.0pt{$\scriptstyle{\mathcal{C}\cup\mathcal{C}_{b}\cup\mathcal{C}_{>}}$\hss}}{\hbox
to0.0pt{$\scriptscriptstyle{\mathcal{C}\cup\mathcal{C}_{b}\cup\mathcal{C}_{>}}$\hss}}}\;\widetilde{G}(r,r_{0},k,m)\mathrm{e}^{-\mathrm{i}kx}dk$
$\displaystyle\;=\;G(x,r;r_{0},m)-I(x)\;=\;\sum_{\mathchoice{\hbox
to0.0pt{\hss$\displaystyle{j:\;\mathrm{Im}(k_{j})<\mathrm{Im}(\mathcal{C})}$\hss}}{\hbox
to0.0pt{\hss$\textstyle{j:\;\mathrm{Im}(k_{j})<\mathrm{Im}(\mathcal{C})}$\hss}}{\hbox
to0.0pt{\hss$\scriptstyle{j:\;\mathrm{Im}(k_{j})<\mathrm{Im}(\mathcal{C})}$\hss}}{\hbox
to0.0pt{\hss$\scriptscriptstyle{j:\;\mathrm{Im}(k_{j})<\mathrm{Im}(\mathcal{C})}$\hss}}}R(k_{j})$
$\displaystyle\text{for}\quad x>0,$ (39b)

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

The poles of $\widetilde{G}$ correspond to zeros of the denominator of
$\widetilde{G}$, as given in (32). 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}\widehat{D}_{2}-\widehat{C}_{2}D_{1}=0$. As discussed
in section 2.5, this occurs when both $\psi_{1}$ and $\psi_{2}$ satisfy both
boundary conditions at $r=0$ and $r=1$. These modal poles can be further
classified into acoustic modes and surface modes: acoustic modes are those for
which $\alpha$ in equation (11) has a small imaginary part, and correspond to
functions which are oscillatory in $r$; and surface modes are those for which
$\alpha$ has a significant imaginary part, and correspond to functions which
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, 2003; Brambley, 2013).

Since the modal poles occur as zeros of
$C_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}=0$, which we shall assume are
simple zeros, the contribution from the residues of these poles are given by

$R(k)=-\mathrm{sgn}(x)\frac{\omega-U(r^{*})k}{2\upi
r^{*}W(r^{*})}\frac{\psi_{1}(\check{r})\psi_{2}(\hat{r})}{\frac{\partial}{\partial
k}\big{(}C_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}\big{)}}\mathrm{e}^{-\mathrm{i}kx}.$
(40)

The second type of poles are the non-modal poles, which occur when
$W(r^{*})=0$. These occur when we loose independence between
$\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ at $r^{*}$. Note from the formula
for $W(r)$, equation (20), that $W(r^{*})=0$ implies that $r^{*}=r_{c}^{+}$.
Since $1-h\leq r^{*}\leq 1$ this can only occur 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 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
(39) above, and only occurs 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 around the critical layer branch cut denoted $\mathcal{C}_{b}$ in figure
3,

$I(x)=\frac{-1}{2\upi}\int_{\mathcal{C}_{b}}\widetilde{G}\mathrm{e}^{-\mathrm{i}kx}\,\mathrm{d}k.$
(41)

However, as it stands, this integral for $I(x)$ is oscillatory, owing to the
$\mathrm{e}^{-\mathrm{i}kx}$ factor in the integrand, and so 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: firstly, it
allows accurate numerical calculation of the integral; secondly, it allows the
derivation of large-$x$ asymptotics using the Method of Steepest Descents; and
thirdly, it brings insight into the various contributions that make up $I(x)$.
In deforming the integration contour, however, we must analytically continue
$\widetilde{G}$ behind the branch cut (as described in section 3.2 above), and
carefully deform around any poles and branch points. This is illustrated
schematically in figure 4.

$\mathrm{Im}(k)$$\mathrm{Im}(k)$$\frac{\omega}{M}$$\frac{\omega}{M}$$k_{<}$$k_{<}$$k_{>}$$k_{>}$$k^{+}$$k^{+}$$k^{-}$$k^{-}$
Figure 4: (Left) 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 $\widetilde{G}^{-}$ and
$\widetilde{G}^{+}$ are denoted $k^{-}$ and $k^{+}$ respectively. (Right) The
integration contour after being transformed onto the steepest decent contour.
Red lines behave as if evaluated below $\frac{\omega}{M}$ (using
$\widetilde{G}^{-}$); blue as if having been analytically continued around the
$\frac{\omega}{M}$ branch point; green as if having been analytically
continued around the $\frac{\omega}{M}$ and $k_{<}$ branch points; and purple
as if analytically continued around all branch points, giving
$\widetilde{G}^{+}$. Note that we have been required to deform contours around
the $k^{+}$ and $k^{-}$ poles.

Note that poles and branch points of $\widetilde{G}$ may exist behind the
critical layer branch cut, and we must therefore use analytic continuations of
$\widetilde{G}$; the reader is reminded that $\widetilde{G}^{+}$ is the
analytic continuation of $\widetilde{G}$ down behind the branch cut from
above, while $\widetilde{G}^{-}$ is the analytic continuation of
$\widetilde{G}$ up behind the branch cut from below. Here, we use the notation
that a pole of $\widetilde{G}^{+}$ with $\mathrm{Re}(k)>\frac{\omega}{M}$ is
denoted $k^{+}$, and a pole of $\widetilde{G}^{-}$ with
$\mathrm{Re}(k)>\frac{\omega}{M}$ is denoted $k^{-}$. Thus, a $k^{+}$ pole
with $\mathrm{Im}(k^{+})<0$ or a $k^{-}$ pole with $\mathrm{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
$\mathrm{Im}(k)<0$, although this is not always the case.

The Steepest Descent contours are where $\mathrm{e}^{-\mathrm{i}kx}$ is
exponentially decaying; 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 $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 $\frac{\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 $\widetilde{p}_{2}(r)$, and the presence of
$\widetilde{p}_{2}(1\\!-\\!h)$, $\widetilde{p}_{2}(r_{0})$ and
$\widetilde{p}_{2}(r)$ in the expression for $\widetilde{G}$; each of these
terms leads to a branch point, respectively at $\frac{\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, $\widetilde{G}$ possesses a pole at $k_{0}$, which is exactly the
non-modal pole referred to above, although there are no poles of
$\widetilde{G}$ at $\frac{\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

$I(x)=I_{\frac{\omega}{M}}(x)+I_{0}(x)+I_{r}(x)+R_{0}^{+}(k_{0})+\sum_{\mathchoice{\hbox
to0.0pt{\hss$\displaystyle{\mathrm{Im}(k^{+})<0}$\hss}}{\hbox
to0.0pt{\hss$\textstyle{\mathrm{Im}(k^{+})<0}$\hss}}{\hbox
to0.0pt{\hss$\scriptstyle{\mathrm{Im}(k^{+})<0}$\hss}}{\hbox
to0.0pt{\hss$\scriptscriptstyle{\mathrm{Im}(k^{+})<0}$\hss}}}R^{+}(k^{+})-\sum_{\mathchoice{\hbox
to0.0pt{\hss$\displaystyle{\mathrm{Im}(k^{-})<0}$\hss}}{\hbox
to0.0pt{\hss$\textstyle{\mathrm{Im}(k^{-})<0}$\hss}}{\hbox
to0.0pt{\hss$\scriptstyle{\mathrm{Im}(k^{-})<0}$\hss}}{\hbox
to0.0pt{\hss$\scriptscriptstyle{\mathrm{Im}(k^{-})<0}$\hss}}}R^{-}(k^{-}),$
(42)

where $R^{\pm}$ is the residue given in (40) evaluated using
$\widetilde{G}^{\pm}$, $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
appendix C.2 and given in (75ck) as

$R_{0}^{+}(k_{0})=\frac{2Mk_{0}^{2}(\omega-
Mk_{0})\mathrm{e}^{-\mathrm{i}k_{0}x}}{3\pi
r_{0}h^{2}\omega(C_{1}^{+}\widehat{D}_{2}-\widehat{C}_{2}D_{1})}\begin{cases}\widehat{D}_{2}\psi_{1}(r)&r<r_{0}\\\
D_{1}\psi_{2}(r)&r>r_{0},\end{cases}$ (43)

the steepest descent integrals are defined as

$I_{q}(x)=\frac{1}{2\upi\mathrm{i}}\int_{0}^{\infty}\Delta\widetilde{G}_{q}(k_{q}-\mathrm{i}\xi)\mathrm{e}^{-\mathrm{i}(k_{q}-\mathrm{i}\xi)x}\,\mathrm{d}\xi,$
(44)

and the jumps across each of the Steepest Descent branch cuts are calculated
in appendix B to be

$\displaystyle\Delta\widetilde{G}_{\frac{\omega}{M}}$
$\displaystyle=\frac{-\left(\omega-U(r^{*})k\right)A}{r^{*}W(r^{*})\big{(}C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}\big{)}}\times\begin{cases}\dfrac{\widehat{D}_{2}^{2}\psi_{1}^{-}(\check{r})\psi_{1}^{-}(\hat{r})}{\left(C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}\right)}&\hat{r}<1\\!-\\!h\\\
\\\
\dfrac{D_{1}\widehat{D}_{2}\psi_{1}^{-}(\check{r})\psi_{2}^{-}(\hat{r})}{\left(C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}\right)}&\check{r}<1\\!-\\!h<\hat{r}\\\
\\\
\dfrac{D_{1}^{2}\psi_{2}^{-}(\check{r})\psi_{2}^{-}(\hat{r})}{\left(C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}\right)}&1\\!-\\!h<\check{r},\end{cases}$
(45a) $\displaystyle\Delta\widetilde{G}_{<}$
$\displaystyle=\frac{-(\omega-U(r^{*})k)}{r^{*}W(r^{*})}\frac{AD_{1}\widetilde{p}_{1}(\check{r})\psi_{2}^{-}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}H\big{(}\check{r}-(1-h)\big{)},$ (45b)
$\displaystyle\Delta\widetilde{G}_{>}$
$\displaystyle=\frac{-(\omega-U(r^{*})k)}{r^{*}W(r^{*})}\frac{A\widehat{D}_{2}\psi_{1}^{-}(\check{r})\widetilde{p}_{1}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}H\big{(}\hat{r}-(1-h)\big{)},$ (45c)

Note that
$\Delta\widetilde{G}_{\frac{\omega}{M}}+\Delta\widetilde{G}_{<}+\Delta\widetilde{G}_{>}=\Delta\widetilde{G}=\widetilde{G}^{+}-\widetilde{G}^{-}$.
While these integrals are now amenable to numerical integration, additional
understanding of the contribution from the three steepest decent 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 (42) contains integrals $I_{q}(x)$
given by (44) which 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, 1918). 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$,

$\frac{1}{2\mathrm{i}\upi}\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\Gamma(\nu+1)\mathrm{e}^{-\mathrm{i}k_{q}x}}{2\mathrm{i}\upi
x^{\nu+1}}+O\big{(}x^{-(\nu+2)}\big{)},$ (46)

where $\Gamma$ is the Gamma function, and in particular, $\Gamma(\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
$\frac{\omega}{k_{q}}$.

In order to find the decay rates of the steepest decent 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 these can be found in appendix C. The
result, given in equations (75ce) and (75cf), is that, for
$k=\frac{\omega}{M}-\mathrm{i}\xi$, as $\xi\to 0$ with $\xi>0$, we find that

$\Delta\widetilde{G}_{\frac{\omega}{M}}\sim\begin{cases}\xi^{3/2}&r_{0}\leq
1-h\\\ \xi^{5/3}&r_{0}>1-h.\end{cases}$ (47)

By Watson’s Lemma, this results in a wave convected with the flow speed
$M=U(1-h)$ and algebraically decaying like $x^{-\frac{5}{2}}$ when the source
is within the region of uniform flow, and $x^{-\frac{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
$\Delta\widetilde{G}_{0}$ as $k\to k_{0}$ is derived in appendix C.2 as

$\Delta\widetilde{G}_{0}=\frac{A\omega
h^{2}U(r_{0})}{6r_{0}Mk_{0}^{2}(r_{0}-1+h)}\frac{(k-k_{0})^{2}}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\upi\mathrm{i}AD_{1}\widehat{D}_{2}}\times\begin{cases}\widehat{D}_{2}\psi_{1}^{-}(r)&r_{0}>r\\\
D_{1}\psi_{2}^{-}(r)&r_{0}<r.\end{cases}$ (48)

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 $\Delta\widetilde{G}_{r}$ as $k\to k_{r}$, it is found in
appendix C.3 that

$\Delta\widetilde{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}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}\times\begin{cases}\widehat{D}_{2}\psi_{1}^{-}(r_{0})&r_{0}<r\\\
D_{1}\psi_{2}^{-}(r_{0})&r_{0}>r.\end{cases}$ (49)

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. (2012a).
We now proceed to compare these results with 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 (1975) for a general flow
profile, and those predicted by Brambley et al. (2012a) for a constant-then-
linear flow profile, in table 1.

| $I_{\frac{\omega}{M}}$ | $I_{r}$ | $I_{0}$ | $R_{0}^{+}(k_{0})$
---|---|---|---|---
| $r_{0}\leq 1-h$ | $r_{0}>1-h$ | $r>1-h$ | $r_{0}>1-h$ | $r_{0}>1-h$
Swinbanks | – | – | – | $x^{-3}$ | 1
Linear BL | $x^{-4}$ | $x^{-5}$ | $x^{-4}$ | $x^{-3}$ | 1
Quadratic BL | $x^{-\frac{5}{2}}$ | $x^{-\frac{7}{2}}$ | $x^{-4}$ | $x^{-3}$ | 1
Table 1: Comparison of the different decay rates given for a general flow
profile by Swinbanks (1975) and for a linear boundary layer flow profile by
Brambley et al. (2012a) against those given here for a quadratic boundary
layer flow profile.

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 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 did not consider the other cases in detail, this work finds
further agreement for the $I_{r}$ contribution with Brambley et al.. In both
the linear and quadratic shear flow cases, when the source is located within
the region of sheared flow, the $I_{0}$ contibution is the slowest decaying
term. When the source is located within the uniform flow region, the
$I_{\frac{\omega}{M}}$ contribution is the slowest decaying term, although
this is matched by $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.

The difference in the behaviour of the $I_{\frac{\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

$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}}.$
(50)

In the case of linear shear flow, the $U^{\prime\prime}$ term is zero for
$k\not=\frac{\omega}{M}$ and the resulting expression is $O(1)$ as
$k\to\frac{\omega}{M}$. In the case of a quadratic shear boundary layer,
$U^{\prime\prime}$ is non zero and dominates $A$ as $k\to\frac{\omega}{M}$,
providing a factor of $(k-\frac{\omega}{M})^{-\frac{1}{2}}$. The remainder of
the differences between the decay rates is explained, for
$I_{\frac{\omega}{M}}$, by the fact that
$(r_{c}^{+}-(1-h))\sim(k-\frac{\omega}{M})^{\frac{1}{2}}$ in the quadratic
shear case, where as for linear shear $(r_{c}-(1-h))\sim(k-\frac{\omega}{M})$.
For the $I_{0}$ and $I_{r}$ contributions, where we do not have $r_{c}^{+}\to
1-h$, and 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 Numerical Results

In this section, the above analysis is illustrated with some numerical
examples. The Frobenius series solutions $\widetilde{p}_{1}$ and
$\widetilde{p}_{2}$ are computed by summing the terms of the series, as given
in appendix A, until a relative error of order $10^{-16}$ is achieved. The
modal poles are found using a variant of the Secant method, and have been
confirmed against results using a finite-difference method applied to the
Pridmore-Brown equation (Brambley, 2011b). When summing the residues of modal
poles, all poles with $|\mathrm{Im}(k)|<400$ have been included.

Throughout this section we show results from four parameter sets, detailed in
table 2.

|  | A1 | A2 | B | C
---|---|---|---|---|---
Frequency |  $\omega$ | 10 | 10 | 31 | 16
Azimuthal order |  $m$ | 5 | 5 | 24 | 24
Centreline Mach number |  $M$ | 0.5 | 0.5 | 0.5 | 0.35
Boundary layer thickness |  $h$ | 0.05 | 0.001 | 0.01 | 0.005
Impedance mass |  $\mu$ | 0.2 | 0.2 | 0.01 | 0.01
Impedance spring |  $K$ | 10 | 10 | 10 | 10
Impedance damper |  $R$ | 2 | 2 | 0.75 | 0.75
Impedance $Z(\omega)$ |  $Z$ |  $2+\mathrm{i}$ |  $2+\mathrm{i}$ |  $0.75-0.01\mathrm{i}$ | $0.75-0.465\mathrm{i}$
Table 2: Parameter sets used for the following numerical results. The
impedance is of mass–spring–damper type,
$Z(\omega)=R+\mathrm{i}\mu\omega-\mathrm{i}K/\omega$.

These parameter sets are inspired by values used in previous studies (Brambley
et al., 2012a; Brambley, 2013; Brambley & Gabard, 2016), and motivated by
application to aeroengine intakes; in particular, parameter set B is intended
to be typical of a rotor-alone tone at takeoff, while parameter set C might
represent the same type of mode during the landing approach.

In section 4.1, we will explore the locations of the modal poles in the
complex $k$-plane. This section will particularly focus on the $k^{\pm}$ modal
poles discussed in 3.3, including tracking these modal poles as they move
behind the critical layer branch cut for certain parameters (by taking
advantage of the ability of the Frobenius series solutions to analytically
continue behind the critical layer branch cut). In section 4.2 we compare the
various contributions from the critical layer branch cut described in section
3.3.3, including their large $x$ behaviour, and show that these agree with the
predicted large $x$ behaviour from section 3.4. In section 4.3 the full
solution in terms of $x$ and $r$ is plotted, and compared these results to the
linear boundary layer case. Finally, in section 4.4, we investigate how the
results vary as we vary the frequency $\omega$, the boundary layer thickness
$h$, and the wall impedance $Z$.

###  4.1 Pole Locations.

The locations of the poles in the complex $k$ plane for parameter sets A and B
are plotted in figure 5.

Figure 5: Location of the poles in the complex $k$ plane for parameter sets A1
(top) and A2 (bottom). Left: For real $\omega$: acoustic modes with
$\mathrm{Re}(k)<\frac{\omega}{M}$ $(\ast)$; $k^{+}$ poles $(+)$; $k^{-}$ poles
$(\times)$; the critical layer branch cut (–); and branch points
$\frac{\omega}{M}$ and $k_{0}$ for $r_{0}=1-\frac{9h}{10}$ $(\circ)$. Right:
Trajectories of poles for $-50<\mathrm{Im}(\omega)<0$. Poles coloured red
(left) and solid lines (right) denote poles contributing to the modal sum.
Poles coloured blue (left) and dashed lines (right) denote poles hidden behind
the branch cut (which varies with $\mathrm{Im}(\omega)$) and do not
contribute.

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 section
3.3.3 above. 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
$\mathrm{Im}(\omega)$ is decreased from zero, following the Briggs-Bers
Criterion (shown in the right-hand plots in figure 5); note that the critical
layer branch cut also moves as a function of $\mathrm{Im}(\omega)$. Of
particular note is that the $k^{+}$ pole for parameter set A1 emerges from
behind the critical layer branch cut as $\mathrm{Im}(\omega)$ is reduced from
zero, becoming a standard modal pole provided $\mathrm{Im}(\omega)$ is taken
sufficiently negative.

As discussed in section 3.3, when the $k^{+}$ 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 it’s contribution
to the modal sum directly, but instead it contributes as part of the branch
cut integral, as seen when deforming onto the steepest decent contour. In this
latter case, we would observe a contribution that decays in $x$. In both
examples, the $k^{-}$ pole in 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 decent 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 section 4.3
below.

Also plotted in figure 5 is the critical layer branch cut for
$k>\frac{\omega}{M}$, together with the non-modal $k_{0}$ pole, which is only
present 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 section.

###  4.2 Branch cut contributions

Figure 6

Figure 6: A comparison of the terms that contribute to the critical layer, for
$r_{0}=1-\frac{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_{\frac{\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).

illustrates the differences between the three types of contributions occurring
due to the presence of the boundary layer: the three steepest descent contour
integrals (figure 6(i)); the $k_{0}$ non-modal pole (figure 6(ii)); and the
$k^{+}$ modal pole when it is located below the branch cut and therefore does
not appear in the modal sum (figure 6(iii)). The non-modal $k_{0}$ pole
appears in these cases to have a small contribution compared to the three
integrals for small $x$, although it is comparable for larger $x$. For small
$x$, the $k^{+}$ poles contribution is far greater than those of the three
integrals and the non-modal $k_{0}$ pole. 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. We can further look at how these
contributions vary as we adjust the location of the source, shown in figure 7.

Figure 7: Plots of the real part of the contribution from integrating around
the branch cut ($\mathrm{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-\frac{9h}{10}$;
(b) $r_{0}=1-\frac{3h}{5}$.

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$. Note that there is no non-modal
$k_{0}$ pole when $r_{0}\leq 1-h$.

Figure 8

Figure 8: Plotted for parameter set A1 are $|I_{\frac{\omega}{M}}|$ (top),
$|I_{r}|$ (middle), and $|I_{0}|$ (bottom). The point source is at
$r_{0}=1-\frac{4h}{5}$ (left) and $r_{0}=0.4$ (right). Lines correspond to
radial locations $r=0.2,0.6,1-\frac{9h}{10}$, and $1-\frac{3h}{5}$. The dashed
line is the predicted far-field rate of decay.

compares the numerically-computed steepest decent integrals with their
predicted far-field rates of decay given in section 3.4, and a good agreement
is seen in all cases.

###  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

(a)(b)

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-\frac{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.

compares a snapshot in the near-field (for small $x$ values) of the wave field
generated by only the stable modal poles (left) with the full solution
including the critical layer and any unstable $k^{+}$ pole (right). 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 plot, 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.

In comparison, figure 10

(a)(b)(c)

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 decent 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-\frac{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$).

shows the behaviour outside the near field of 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 11

Figure 11: Plotting the real values of the full solution for a quadratic
boundary layer flow profile (top) and a linear boundary layer flow profile
(bottom) (from Brambley et al., 2012a). From left to right, parameters are:
set A1 with $r_{0}=0.4$; set A1 with $r_{0}=1-\frac{4h}{5}=0.96$; set A2 with
$r_{0}=0.4$; and set A2 with $r_{0}=1-\frac{4h}{5}=0.9992$.

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.,
2012a). 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 only found to be unstable 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 rate of the
instability can be significantly different.

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 next.

###  4.4 Variation of results with changing parameters

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

Figure 12

Figure 12: Motion of the modal poles for parameter set C as one parameters is
varied (arrows show the motion as the parameter is increased): (a) varying $h$
in $(0.001,0.5)$; (b) varying $\mathrm{Im}(Z)$ in $(-\infty,\infty)$, with a
dot showing hard-walled values; (c) varying $\mathrm{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$ are varied.

illustrates how the $k^{+}$ and $k^{-}$ modal poles vary with boundary layer
thickness, frequency, impedence 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 $\mathrm{Im}(Z)$ being unstable and
nearly hard-walled values of $|\mathrm{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 either to the modal
sum or the critical layer branch cut.

##  5 Conclusions

In this work we have considered a cylindrical duct containing a parallel mean
flow that is uniform everywhere except within a boundary layer of thickness
$h$ near the wall. Within the boundary layer, which need not be thin, the flow
has a quadratic profile and satisfies the non-slip boundary condition at the
duct wall, whilst maintaining a $C^{1}$ continuous flow. For such a flow
profile for $0<h<1$, irrespective of the the thickness of the boundary layer,
a solution to the Pridmore-Brown equation has been constructed making use of
two Frobenius series expansions, valid for any wave number $k$. This enables
the evaluation of the Greens function of the Pridmore-Brown equation,
equivalent to the solution for a point mass source in the linearised Euler
equations. The Green’s function solution is found to consist of a sum of the
usual acoustic duct modes, plus a non-modal contribution from the critical
layer branch cut. Full source code is provided in the supplementary material
to evaluate all solutions presented here.

The Frobenius series method employed here has two particular advantages over
other numerical methods to solve the Pridmore-Brown equation (such as finite
differences, e.g., Brambley, 2011b). The first is that the Frobenius series,
being a series solution about the critical points of the equation, is at its
most accurate near the critical layer singularity found in the Pridmore-Brown
equation. This allows for accurate numerical solutions near the critical
layer, required for the integration around the critical layer singularity and
its associated branch cut to evaluate their effect on the resulting pressure
field. Other numerical methods such as finite difference are typically at
their least accurate near the critical layer (Brambley et al., 2012a).
Moreover, the Frobenius series solution makes explicit the branch cut along
the critical layer, allowing for analytic continuation of the solution behind
the branch cut. This allows for tracking hydrodynamic instability surface wave
modes as they become stable and enter the critical layer (as seen figure 12),
which makes it significantly easier to track the boundary between stable and
unstable behaviour.

An advantage of considering this particular quadratic flow profile is that the
origins of the critical layer are on a more solid footing. For the linear flow
profile (Brambley et al., 2012a), the critical layer was due to the
cylindrical geometry of the duct, where as in general the critical layer is
due to a non-zero second derivative of the sheared flow profile. This also
allows comparison to previous works, such as that of Swinbanks (1975) and
(Félix & Pagneux, 2007). Further, as the quadratic flow profile has a
continuous first derivative, we have also been able to investigate the
specific case of a point mass source at the boundary between uniform and
sheared flow, $r_{0}=1-h$, and we find this case retains the same behaviour as
a point mass source within the region of uniform flow. In contrast, for the
linear flow profile, $r_{0}\to 1-h$ is a singular limit. We therefore believe
the results of the quadratic boundary layer flow profile to be in some ways
typical of solutions for a general boundary layer profile.

The final solution for the Green’s function for a point mass source was found
to consist of a number of contributions. This solution is dominated, both
upstream and, in the near field, downstream too, by the sum of modal poles.
The modal poles, including acoustic and surface modes, are well known, and are
typically used in mode-matching numerical methods. One complication found here
to the surface modes is that a particular surface mode, here labelled $k^{+}$,
is found to sometimes disappear behind the critical layer branch cut (or, in
other words, into the continuous spectrum). The contribution of this mode is
not lost, however, and is in effect added to the critical layer contribution.
In general, the modal poles present difficulty only in establishing which
poles contribute upstream $(x<0)$ or downstream $(x>0)$ of the source. This
can be established through application of the Briggs-Bers criterion, as
summarised in section 3.3.

The effect of the critical layer, the focus of this work, always contributes
downstream of the source, and is the dominant contribution to the far-field
pressure downstream of the point mass source. This contribution, which results
from integrating the Fourier inversion contour around the critical layer
branch cut, may be viewed in three parts. The first contribution is from the
$k_{0}$ non-modal pole, only present when $r_{0}>1-h$, which does not decay
with distance from the point source and is therefore dominant in the farfield
downstream of the source. This contribution is similar to that described in
the linear flow profile case (Brambley et al., 2012a), and may be interpreted
similarly as a hydrodynamic vorticity wave generated from the point mass
source interacting with the sheared mean flow, traveling downstream with phase
velocity equal to the local fluid velocity $U(r_{0})$. The second contribution
to the critical layer is from the steepest-descent non-oscillatory integrals
$I_{\frac{\omega}{M}}$, $I_{r}$ and $I_{0}$, which are the results of
accounting for the branch points coming from the critical points of the
Pridmore-Brown equation. These contributions have a phase speed equal to,
respectively, the uniform flow speed, $U(r)$, and $U(r_{0})$, and decay
algebraically in the far-field downstream of the point source. The final
contribution is from any modal pole that is hiding behind the branch cut, such
as from a $k^{+}$ surface wave mode that has stabilized by moving into the
critical layer branch cut from above. These poles, while looking very much
like ordinary modal poles, are not able to be found by traditional numerical
methods, as they require analytically continuing behind the critical layer
branch cut. While these poles decay exponentially with distance downstream of
the point source, their decay rate may be slower than any other acoustic modal
pole, depending on the parameters used, and so may still be significant in the
far-field; this was seen for parameter set C in figure 10.

The $k^{+}$ modal pole may be present as a hydrodynamic instability surface
wave, or as a stable mode included within the critical layer branch cut
contribution. Interestingly, in the linear flow profile case (Brambley et al.,
2012a), this mode was always an instability and was never hidden behind the
critical layer branch cut. From the results of figure 12, we expect that this
mode is stable for quadratic flow boundary layer profiles when the boundary
layer is sufficiently thick or the frequency is sufficiently high, although
the specific stability boundary also depends on the impedance $Z$.

For the linear flow profile boundary layer (Brambley et al., 2012a), a $k^{-}$
pole was found below and behind the critical layer branch cut that contributed
to the critical layer. For the quadratic flow profile boundary layer here,
this $k^{-}$ pole is always found to be above the critical layer branch cut,
and so never contributes. We believe that this $k^{-}$ pole was an artifact of
the unphysical linear boundary layer profile, although we have no direct way
of demonstrating this. Incidentally, for the linear flow profile boundary
layer, Brambley et al. (2012a) argued that the $k^{+}$ pole could never be
behind the critical layer branch cut, as this would cause an unphysical
discontinuity in the final solution; in fact, it is found here that when the
$k^{+}$ pole is behind the branch cut, the unphysical discontinuity in the
$k^{+}$ pole contribution is exactly cancelled by the $I_{r}$ steepest descent
contour contribution, resulting in a continuous solution, as expected.

The various decay rates of the components of the critical layer have
previously been predicted by Swinbanks (1975) and Brambley et al. (2012a); and
a summary can be found in table 1. Swinbanks (1975) only considered the
contribution from waves with phase velocity $U(r_{0})$, which are only present
for a point mass source within the region of sheared flow, $r_{0}>1-h$.
Swinbanks (1975) predicted these to behave as a constant amplitude plus a
decay as $O(x^{-3})$ in the far field. Brambley et al. (2012a) found the same
result, despite Swinbanks considering a two dimensional flow in a rectilinear
duct with an arbitrary flow profile and Brambley et al. considering only a
constant-then-linear flow profile in a three dimensional cylindrical duct; in
particular Swinbanks emphasised the importance of the non-zero second
derivative of the mean flow profile, which is identically zero for a constant-
then-linear flow profile. As a result, it would not have been unsurprising if
these results were different. Here, the same result is again found, with the
constant amplitude coming from the $k_{0}$ non-modal pole and the algebraic
decay coming from the $I_{0}$ integral. This shows that this agreement is not
by coincidence. For the critical layer contribution that propagates with phase
velocity $U(r)$ when $r$ is within the boundary layer, we also find here the
same result given by Brambley et al. (2012a) of an $O(x^{-4})$ far-field
decay.

The critical layer also contributes a term with phase velocity equal to the
uniform flow velocity $M$, which is always present, and which dominates the
other critical layer contributions in the far field whenever the point mass
source is in the uniform flow region, $r_{0}<1-h$. The amplitude of this term
can decay at two different rates depending on the location of the source. When
the source is within the uniform flow a decay rate of $O(x^{-\frac{5}{2}})$ is
found, while when the source is within the sheared flow we instead have a
faster rate of decay of $O(x^{-\frac{7}{2}})$. These results differ from those
found by Brambley et al. (2012a) in the linear flow profile boundary layer
case, despite corresponding to the same physical behaviour. This difference
may be understood as result of both the difference in the overall shape of the
flow profile, and the importance of the second derivative of the flow. Indeed,
we conjecture that these scalings will differ depending on the flow profile
within the boundary layer, and an example discussion of this for
$n$-polynomial flow profiles is given in appendix D.

In most aeroacoustic analyses, particularly those involving mode matching, the
critical layer is either implicitly or explicitly ignored. The work here
suggests that this may be valid in the near-field provided not all acoustic
modes are cut-off, although even in the near-field the critical layer can be
dominant if all acoustic modes are cut-off, as shown in figure 9 for parameter
set C. However, it is certainly not valid to ignore the critical layer
downstream in the far-field, when the critical layer will be the dominant
contribution. Moreover, without considering the critical layer, it would not
be apparent whether an barely-stable hydrodynamic surface wave is present only
just hidden behind the critical layer branch cut (or, in other words, within
the continuous spectrum).

There are a number of possible avenues for further investigation following on
from this work. One of practical importance concerns whether the hydrodynamic
surface wave $k^{+}$ can be accurately predicted using a surface wave
dispersion relation (e.g. Brambley, 2013), especially when the $k^{+}$ pole is
located behind the critical layer branch cut; our experience in this work has
been that it cannot, at least with the simplified surface wave dispersion
relations that assume a thin boundary layer with a linear flow profile,
although more complicated surface wave dispersion relations may prove more
accurate. Another possibility for further investigation is to consider
parameters on the stability boundary when the $k^{+}$ pole is neutrally
stable. In this case, the $k^{+}$ pole is exactly located on the critical
layer branch cut, and there would also exist a value of $r_{0}$ for which the
non-modal $k_{0}$ pole and the $k^{+}$ pole overlap; this case has been
explicitly excluded here. While this may seem a rather contrived case, a
distributed sound source would correspond to an integral of source strength
over all values of $r_{0}$, and so $k_{0}$ and $k^{+}$ coinciding could be
expected to occur for any parameters for leading to exact neutral stability.
Finally, the critical layer may be regularized by considering either viscosity
or weak nonlinearity, and it would be interesting to investigate how the
results presented here are recovered in the inviscid or small-amplitude
limits. In particular, for viscous thin boundary layers, the critical layer is
recovered as a caustic in the high-frequency limit (Brambley, 2011a).

###### Acknowledgements.

MJK was supported in this work through the University of Warwick MASDOC
Doctoral Training Centre, and gratefully acknowledges their support. EJB
gratefully acknowledges the support of a Royal Society University Research
Fellowship (UF150695 and RGF\EA\180284). RL as supported in this work through
a research internship funded by the Royal Society (RGF\EA\180284), and would
also like to thank the CAPES Foundation, Ministry of Education of Brazil for
the award of a BRAFITEC scholarship. The contribution of MR was supported by
an EPSRC UROP undergraduate research summer internship (2015, DAMTP,
University of Cambridge).

##  Appendix A Frobenius series solutions to the Pridmore-Brown equation witha
quadratic mean flow profile

In this appendix, we use a Frobenius expansion method to solve the homogeneous
Pridmore-Brown equation (9),

$\widetilde{p}^{\prime\prime}+\left(\frac{2kU^{\prime}}{\omega-U(r)k}+\frac{1}{r}\right)\widetilde{p}^{\prime}+\left((\omega-U(r)k)^{2}-k^{2}-\frac{m^{2}}{r^{2}}\right)\widetilde{p}=0,$
(51)

for the flow profile (8),

$U(r)=\begin{cases}M&0\leq r\leq 1-h\\\ M(1-(1-\frac{1-r}{h})^{2})&1-h\leq
r\leq 1,\end{cases}$ (52)

in the quadratic flow region $1-h\leq r\leq 1$. The Pridmore-Brown equation
(51) has regular singularities at $r=0$ and at $r=r_{c}$, where
$\omega-U(r_{c})k=0$. For the quadratic flow profile (52), the solutions of
$\omega-U(r_{c})k=0$ are given by (12),

$\displaystyle r_{c}^{\pm}$ $\displaystyle=1-h\pm Q,$
$\displaystyle\text{where}\quad Q$
$\displaystyle=h\sqrt{1-\frac{\omega}{Mk}}.$ (53)

This results in the Pridmore-Brown equation in the quadratic flow region
$1-h\leq r\leq 1$ being given by

$\widetilde{p}^{\prime\prime}+\left(\frac{2}{r-r_{c}^{+}}+\frac{2}{r-r_{c}^{-}}+\frac{1}{r}\right)\widetilde{p}^{\prime}+\left(\frac{M^{2}k^{2}}{h^{4}}(r-r_{c}^{+})^{2}(r-r_{c}^{-})^{2}-k^{2}-\frac{m^{2}}{r^{2}}\right)\widetilde{p}=0.$
(54)

We choose $\mathrm{Re}(Q)\geq 0$ and consider the Frobenius expansion about
$r=r_{c}^{+}$.

###  A.1 Frobenius expansion about $r=r_{c}^{+}$

Following Brambley et al. (2012a), we propose a Frobenius expansion about the
regular singularity $r_{c}^{+}$,

$\displaystyle\widetilde{p}(r)$
$\displaystyle=\sum_{n=0}^{\infty}a_{n}(r-r_{c}^{+})^{n+\sigma}$
$\displaystyle\text{with}\qquad a_{0}\neq 0.$ (55)

We substitute (55) into (54) and expand using a Laurant series. Specifying
that $a_{0}\neq 0$ results in the requirement that $\sigma(\sigma-3)=0$. By
Fuchs theorem (Teschl, 2012), this gives a pair of linearly independent
solutions of the form

$\displaystyle\widetilde{p}_{c1}(r)$
$\displaystyle=\sum_{n=0}^{\infty}a_{n}(r-r_{c}^{+})^{n+3},$ (56a)
$\displaystyle\widetilde{p}_{c2}(r)$
$\displaystyle=A\widetilde{p}_{c1}(r)\log(r-r_{c}^{+})+\sum_{n=0}^{\infty}b_{n}(r-r_{c}^{+})^{n}.$
(56b) The coefficients $a_{n}$ and $b_{n}$ are then given by setting the
remaining terms of the Laurent expansion of (54) to zero, resulting in the
recurrence relation $\displaystyle a_{n}$
$\displaystyle=\frac{1}{n(n+3)}\Bigg{[}k^{2}a_{n-2}-\frac{k^{2}M^{2}}{h^{4}}\left(a_{n-6}+4Qa_{n-5}+4Q^{2}a_{n-4}\right)$
$\displaystyle\qquad\qquad\qquad-\sum_{j=0}^{n-1}(-1)^{j}(n+2-j)a_{n-1-j}\left(\frac{1}{(r_{c}^{+})^{j+1}}-\frac{2}{(2Q)^{j+1}}\right)$
$\displaystyle\qquad\qquad\qquad+m^{2}\sum_{j=0}^{n-2}\frac{(-1)^{j}}{(r_{c}^{+})^{j+2}}(j+1)a_{n-2-j}\Bigg{]},$
(56c) $\displaystyle b_{n}$
$\displaystyle=-\frac{1}{n(n-3)}\Bigg{[}A\Bigg{(}\\!(2n-3)a_{n-3}+\sum_{j=0}^{n-4}a_{n-4-j}(-1)^{j}\\!\left(\frac{1}{(r_{c}^{+})^{j+1}}-\frac{2}{(2Q)^{j+1}}\right)\\!\\!\Bigg{)}$
$\displaystyle\qquad\qquad\qquad-k^{2}b_{n-2}+\frac{k^{2}M^{2}}{h^{4}}\left(b_{n-6}+4Qb_{n-5}+4Q^{2}b_{n-4}\right)$
$\displaystyle\qquad\qquad\qquad+\sum_{j=0}^{n-1}(-1)^{j}(n-1-j)b_{n-1-j}\\!\left(\frac{1}{(r_{c}^{+})^{j+1}}-\frac{2}{(2Q)^{j+1}}\right)$
$\displaystyle\qquad\qquad\qquad-m^{2}\sum_{j=0}^{n-2}\frac{(-1)^{j}}{(r_{c}^{+})^{j+2}}(j+1)b_{n-2-j}\Bigg{]},$
(56d)

where we take $a_{n}=b_{n}=0$ for $n<0$. Requiring $a_{0}=b_{0}=1$, we then
find that

$\displaystyle b_{1}$ $\displaystyle=0,$ $\displaystyle b_{2}=$
$\displaystyle-\frac{1}{2}\\!\left(\\!k^{2}+\left(\frac{m}{r_{c}^{+}}\right)^{\\!\\!2}\right)\\!,$
(57)

and that $b_{3}$ is arbitrary, so we choose $b_{3}=0$. However, for the
recurrence relation involving $b_{3}$ on the left hand side to hold, we also
require that $A$ is chosen to be

$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}}.$
(58)

Here, the notation $\widetilde{p}_{c1}$ and $\widetilde{p}_{c2}$ denotes that
these are two linearly independent solutions for $\widetilde{p}$ about the
critical point $r_{c}^{+}$.

The Frobenius series solutions (56) are limited by a radius of convergence, in
that the series converge if $|r-r_{c}^{+}|<R$ for some radius of convergence
$R$. This $R$ is the distance between $r_{c}^{+}$ and the next nearest
singularity of the Pridmore-Brown equation, which is either at $r=0$ or at
$r=r_{c}^{-}$, and hence

$R=\min\big{\\{}|1-h+Q|,2|Q|\big{\\}}.$ (59)

The choice of $r_{c}^{+}$ as the singularity to expand around means that this
expansion maximizes the region of $[1-h,1]$ contained within the radius of
convergence. This is shown schematically in figure 13.

$0$$\mathrm{Re}(r)$$\mathrm{Im}(r)$$1-h$$1$$r_{c}^{+}$$r_{c}^{-}$(a)$0$$\mathrm{Re}(r)$$\mathrm{Im}(r)$$1$$r_{c}^{+}$$r_{c}^{-}$(b)
Figure 13: Schemetic of possible locations of the $r_{c}^{\pm}$ critical
points in the complex $r$-plane. (a) The radius of convergence of the
expansion about $r_{c}^{+}$ covers the region of interest $r\in[1-h,1]$. (b)
The radius of convergence about $r_{c}^{+}$ is insufficient to cover
$r\in[1-h,1]$ (as $k$ is too close to $\frac{\omega}{M}$).

It can be observed that these solutions are not always valid for all of
$r\in[1-h,1]$. In particular, in the case of $k\to\frac{\omega}{M}$ we observe
$r_{c}^{\pm}\to(1-h)$ and the radius of convergence $R\to 0$.

###  A.2 Frobenius expansion about $r=1$

In order to cover the remainder of the domain $[1-h,1]$, we construct a second
series solution about $r=1$,

$\displaystyle\widetilde{p}_{11}(r)$
$\displaystyle=\sum_{n=0}^{\infty}\alpha_{n}(r-1)^{(n+1)},$
$\displaystyle\widetilde{p}_{12}(r)$
$\displaystyle=\sum_{n=0}^{\infty}\beta_{n}(r-1)^{n}.$ (60)

Specifying that $\alpha_{0}=\beta_{0}=1$, this results in the recurrence
relation

$\displaystyle\alpha_{n}=-\frac{1}{n(n+1)}\Bigg{[}$
$\displaystyle-k^{2}\alpha_{n-2}+\frac{k^{2}M^{2}}{h^{4}}\bigg{(}\alpha_{n-6}+4h\alpha_{n-5}+2(3h^{2}-Q^{2})\alpha_{n-4}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+4h(h^{2}-Q^{2})\alpha_{n-3}+(h^{2}-Q^{2})^{2}\alpha_{n-2}\bigg{)}$
$\displaystyle+\sum_{j=0}^{n-1}(-1)^{j}(n-j)\alpha_{n-j-1}\left(1-\frac{2}{(h+Q)^{j+1}}-\frac{2}{(h-Q)^{j+1}}\right)$
$\displaystyle-m^{2}\sum_{j=0}^{n-2}(-1)^{j}(j+1)\alpha_{n-j-2}\Bigg{]}.$
(61a) $\displaystyle\beta_{n}=-\frac{1}{n(n-1)}\Bigg{[}$
$\displaystyle-k^{2}\beta_{n-2}+\frac{k^{2}M^{2}}{h^{4}}\bigg{(}\beta_{n-6}+4h\beta_{n-5}+2(3h^{2}-Q^{2})\beta_{n-4}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+4h(h^{2}-Q^{2})\beta_{n-3}+(h^{2}-Q^{2})^{2}\beta_{n-2}\bigg{)}$
$\displaystyle+\sum_{j=0}^{n-1}(-1)^{j}(n-j-1)\beta_{n-j-1}\left(1-\frac{2}{(h+Q)^{j+1}}-\frac{2}{(h-Q)^{j+1}}\right)$
$\displaystyle-m^{2}\sum_{j=0}^{n-2}(-1)^{j}(j+1)\beta_{n-j-2}\Bigg{]}.$ (61b)

with $\alpha_{n}=\beta_{n}=0$ for $n<0$. These solutions are labelled
$\widetilde{p}_{11}$ and $\widetilde{p}_{12}$ to indicate they are two
linearly independent solutions to $\widetilde{p}$ expanded about the point
$r=1$.

###  A.3 A homogeneous solution valid across $[1-h,1]$

We now construct solutions to the homogeneous Pridmore-Brown equation
$\widetilde{p}_{1}(r)$ and $\widetilde{p}_{2}(r)$ such that they are valid
across the whole of $[1-h,1]$. We set

$\displaystyle\widetilde{p}_{1}(r)$
$\displaystyle=\begin{cases}\widetilde{p}_{c1}(r)&|r-r_{c}^{+}|<R\\\
A_{1}\widetilde{p}_{11}(r)+B_{1}\widetilde{p}_{12}(r)&\text{otherwise}\end{cases}$
(62a) $\displaystyle\widetilde{p}_{2}(r)$
$\displaystyle=\begin{cases}\widetilde{p}_{c2}(r)&|r-r_{c}^{+}|<R\\\
A_{2}\widetilde{p}_{11}(r)+B_{2}\widetilde{p}{12}(r)&\text{otherwise}\end{cases}$
(62b)

First of all, note that these expansions are sufficient for a uniformly-valid
expansion, as sketched in figure 14.

$0$$\mathrm{Re}(r)$$\mathrm{Im}(r)$$1$$r_{c}^{+}$$r_{c}^{-}$ Figure 14: As for
figure 13(b) with the radius of convergence for $\widetilde{p}_{11}$ and
$\widetilde{p}_{12}$ added.

Note also from figure 14 that the regions of convergence of the
$\widetilde{p}_{c}$ solutions and the $\widetilde{p}_{1}$ solutions always
overlap (except when $k=\frac{\omega}{M}$, which we exclude here). For any
real $\bar{r}>\mathrm{Re}(r_{c}^{+})$ contained within both regions of
convergence, we may find the coefficients $A_{1}$, $A_{2}$, $B_{1}$ and
$B_{2}$ are found by requiring continuity and continuous derivatives at
$r=\bar{r}$:

$\left[\begin{array}[]{cc}A_{1}&A_{2}\\\
B_{1}&B_{2}\end{array}\right]=\left[\begin{array}[]{cc}\widetilde{p}_{11}(\bar{r})&\widetilde{p}_{12}(\bar{r})\\\
\widetilde{p}_{11}^{\prime}(\bar{r})&\widetilde{p}_{12}^{\prime}(\bar{r})\end{array}\right]^{-1}\left[\begin{array}[]{cc}\widetilde{p}_{c1}(\bar{r})&\widetilde{p}_{c2}(\bar{r})\\\
\widetilde{p}_{c1}^{\prime}(\bar{r})&\widetilde{p}_{c2}^{\prime}(\bar{r})\end{array}\right]$
(63)

These coefficients $A_{1},B_{1},A_{2}$ and $B_{2}$ are independent of the
specific choice of $\bar{r}$, and the resulting solutions $\widetilde{p}_{1}$
and $\widetilde{p}_{2}$ have not only $C^{1}$ continuity but $C^{\infty}$,
since both are solutions to the Pridmore-Brown equation. In effect,
$\widetilde{p}_{1}$ analytically continues $\widetilde{p}_{c1}$ beyond its
radius of convergence, and similarly $\widetilde{p}_{2}$ analytically
continues $\widetilde{p}_{c2}$.

As described in (17), there is a jump in $\widetilde{p}_{c2}$ across the
critical layer branch cut due to the $\log$ term in (56b). If the radius of
convergence $R$ is sufficiently large that $r=1$ is within the radius of
convergence, then no matching coefficients are needed, and this jump in
$\widetilde{p}_{c2}$ obviously carries through to $\widetilde{p}_{2}$. In the
other case, that $R$ is sufficiently small that matching is needed, it follows
that $\bar{r}<1$. In this case, there is no jump in the matching coefficients
$A_{1}$, $A_{2}$, $B_{1}$ and $B_{2}$ as $\bar{r}>\mathrm{Re}(r_{c}^{+})$, and
hence

$\Delta\widetilde{p}_{2}(r)=-2\upi\mathrm{i}A\widetilde{p}_{1}(r)H(r_{c}^{+}-r).$
(64)

This is analogous to the jump in $\widetilde{p}_{2}$ given in (17), and shows
that the jump in $\widetilde{p}_{c2}$ carries through the analytic
continuation, as might have been expected a priori.

###  A.4 The Wronskian of $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$

We define the Wronskian of $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ to be

$W(r)=\mathcal{W}(\widetilde{p}_{1},\widetilde{p}_{2};r)=\widetilde{p}_{1}(r)\widetilde{p}_{2}^{\prime}(r)-\widetilde{p}_{2}(r)\widetilde{p}_{1}^{\prime}(r).$
(65)

Since $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ are solutions to the
homogeneous Pridmore-Brown equation (51), we have that

$\displaystyle W^{\prime}+\\!\left(\frac{2kU^{\prime}}{\omega-
Uk}+\frac{1}{r}\right)\\!W$ $\displaystyle=0$ $\displaystyle\Rightarrow$
$\displaystyle W(r)$
$\displaystyle\propto\frac{(r-r_{c}^{+})^{2}(r-r_{c}^{-})^{2}}{r}.$ (66)

By considering the behaviour of $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ as
$r\to r_{c}^{+}$, we find that
$W(r)=-3(r-r_{c}^{+})^{2}+O\big{(}(r-r_{c}^{+})^{4}\big{)}$, and so the
constant of proportionality can be found, giving

$W(r)=-\frac{3}{4}\frac{r_{c}^{+}(r-r_{c}^{+})^{2}(r-r_{c}^{-})^{2}}{rQ^{2}}.$
(67)

##  Appendix B The jump in $\widetilde{G}$ across the critical layer branch
cut

In this appendix, we split the jump in $\widetilde{G}$ across the critical
layer branch cut, $\Delta\widetilde{G}$, into its various components about the
three possible branch points $\frac{\omega}{M}$, $k_{0}$ and $k_{r}$. For this
reason, we restrict attention to $k\in[\frac{\omega}{M},\infty)$, that is, to
$k$ on the critical layer branch cut. In this case,
$r_{c}^{+}(k)\in[1\\!-\\!h,1)$, and $r_{c}^{+}(k)$ is an increasing function
of $k$. Recall from (34) that

$\displaystyle\Delta\widetilde{G}=$
$\displaystyle-\frac{\omega-U(r^{*})k}{2\mathrm{i}\pi
r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}$ (68)
$\displaystyle\quad\times\left[\frac{2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}\psi^{-}_{1}(\check{r})\psi^{-}_{2}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}-\psi^{-}_{1}(\check{r})\Delta\psi_{2}(\hat{r})-\Delta\psi_{1}(\check{r})\psi^{-}_{2}(\hat{r})-\Delta\psi_{1}(\check{r})\Delta\psi_{2}(\hat{r})\right],$

with $\Delta\psi_{1}$ and $\Delta\psi_{2}$ given in (25) as

$\displaystyle\Delta\psi_{1}(r)$ $\displaystyle=\begin{cases}0&r<r_{c}^{+}\\\
2\mathrm{i}\upi AD_{1}\widetilde{p}_{1}&r\geq r_{c}^{+},\end{cases}$ (69a)
$\displaystyle\Delta\psi_{2}(r)$
$\displaystyle=\begin{cases}\Delta\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2}H_{m}^{(1)}(\alpha
r)+\Delta\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2}H_{m}^{(2)}(\alpha
r)&0\leq r\leq 1-h\\\ -2\mathrm{i}\upi
A\widetilde{p}_{1}(r)\widehat{D}_{2}&1-h\leq r\leq r_{c}^{+}\\\
0&r_{c}^{+}<r\leq 1.\end{cases}$ (69b)

Note that since $\check{r}<\hat{r}$ it must be that
$\Delta\psi_{1}(\check{r})\Delta\psi_{2}(\hat{r})=0$ in all cases.

When $r,r_{0}<1-h$ then for any $k$ on the branch cut we have that
$\Delta\psi_{1}=0$ and $\Delta\psi_{2}\neq 0$. This means that we have the
same formula for $\Delta\widetilde{G}$ for any $k$ on the branch cut in this
case, so that $\frac{\omega}{M}$ is the only branch point of
$\Delta\widetilde{G}$. Hence, we write
$\Delta\widetilde{G}=\Delta\widetilde{G}_{\frac{\omega}{M}}$, where

$\displaystyle\Delta\widetilde{G}_{\frac{\omega}{M}}=$
$\displaystyle-\frac{\omega-U(r^{*})k}{2\mathrm{i}\pi
r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}$ $\displaystyle\quad\times\left[\frac{2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}\psi^{-}_{1}(\check{r})\psi^{-}_{2}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}-\psi^{-}_{1}(\check{r})\Delta\psi_{2}(\hat{r})\right],$
$\displaystyle=$
$\displaystyle-\frac{\omega-U(r^{*})k}{r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}\frac{A\widehat{D}_{2}^{2}\psi^{-}_{1}(\check{r})\psi^{-}_{1}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}$
(70a)

When $\check{r}<1-h<\hat{r}$ then the formula for $\Delta\widetilde{G}$
depends on whether $\frac{\omega}{M}<k<k_{>}$ or $k>k_{>}$. In this case, we
set $\Delta\widetilde{G}=\Delta\widetilde{G}_{\frac{\omega}{M}}$ for
$\frac{\omega}{M}_{k}<k_{>}$, and
$\Delta\widetilde{G}=\Delta\widetilde{G}_{\frac{\omega}{M}}+\Delta\widetilde{G}_{>}$
for $k>k_{>}$, so that $\Delta\widetilde{G}_{>}$ is the correction required
for $k>k_{>}$. In effect, $\Delta\widetilde{G}$ has two branch points, one at
$\frac{\omega}{M}$ and one at $k_{>}$ in this case, and by making this
definition we may write

$\int_{\frac{\omega}{M}}^{\infty}\Delta\widetilde{G}\mathrm{e}^{-\mathrm{i}kx}\,\mathrm{d}k=\int_{\frac{\omega}{M}}^{\infty}\Delta\widetilde{G}_{\frac{\omega}{M}}\mathrm{e}^{-\mathrm{i}kx}\,\mathrm{d}k+\int_{k_{>}}^{\infty}\Delta\widetilde{G}_{>}\mathrm{e}^{-\mathrm{i}kx}\,\mathrm{d}k.$
(70b)

By considering (68) in this case, we find that

$\displaystyle\Delta\widetilde{G}_{\frac{\omega}{M}}=$
$\displaystyle-\frac{\omega-U(r^{*})k}{r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}\frac{AD_{1}\widehat{D}_{2}\psi^{-}_{1}(\check{r})\psi^{-}_{2}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}},$
(70c) $\displaystyle\Delta\widetilde{G}_{>}=$
$\displaystyle\phantom{{}-}\frac{\omega-U(r^{*})k}{2\mathrm{i}\pi
r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}\psi^{-}_{1}(\check{r})\Delta\psi_{2}(\hat{r}),$
$\displaystyle=$
$\displaystyle-\frac{\omega-U(r^{*})k}{r^{*}W(r^{*})}\frac{A\widehat{D}_{2}\psi^{-}_{1}(\check{r})\widetilde{p}_{1}(\hat{r})}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}$ (70d)

Finally, when we have $1-h<\check{r}$ we must consider three cases:
$\frac{\omega}{M}<k<k_{<}$, $k_{<}<k<k_{>}$ and $k_{>}<k$. Similarly to the
previous case, we consider
$\Delta\widetilde{G}_{\frac{\omega}{M}}=\Delta\widetilde{G}$ for
$\frac{\omega}{M}<k<k_{<}$, and take $\Delta\widetilde{G}_{<}$ and
$\Delta\widetilde{G}_{>}$ to be correction terms as $k$ crosses $k_{<}$ and
$k_{>}$ respectively. This leads to

$\displaystyle\Delta\widetilde{G}_{\frac{\omega}{M}}=$
$\displaystyle-\frac{\omega-U(r^{*})k}{2\mathrm{i}\pi
r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}$ $\displaystyle\quad\times\left[\frac{2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}\psi^{-}_{1}(\check{r})\psi^{-}_{2}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}-\Delta\psi^{-}_{1}(\check{r})\psi_{2}(\hat{r})\right],$
$\displaystyle=$
$\displaystyle-\frac{\omega-U(r^{*})k}{r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}\frac{AD_{1}^{2}\psi^{-}_{2}(\check{r})\psi^{-}_{2}(\hat{r})}{C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}$
(70e) $\displaystyle\Delta\widetilde{G}_{<}=$
$\displaystyle-\frac{\omega-U(r^{*})k}{2\mathrm{i}\pi
r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}\Delta\psi^{-}_{1}(\check{r})\psi_{2}(\hat{r}),$
$\displaystyle=$
$\displaystyle-\frac{\omega-U(r^{*})k}{r^{*}W(r^{*})}\frac{AD_{1}\widetilde{p}_{1}(\hat{r})\psi^{-}_{2}(\hat{r})}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}$ (70f) $\displaystyle\Delta\widetilde{G}_{>}=$
$\displaystyle\phantom{-}\frac{\omega-U(r^{*})k}{2\mathrm{i}\pi
r^{*}W(r^{*})}\frac{1}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}\psi^{-}_{1}(\check{r})\Delta\psi_{2}(\hat{r}),$
$\displaystyle=$
$\displaystyle-\frac{\omega-U(r^{*})k}{r^{*}W(r^{*})}\frac{A\widehat{D}_{2}\psi^{-}_{1}(\check{r})\widetilde{p}_{1}(\hat{r})}{C^{-}_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}}$ (70g)

##  Appendix C Asymptotic behaviours of $\widetilde{G}$ and
$\Delta\widetilde{G}_{q}$

In order to find the residue contribution of the non-modal pole at $k=k_{0}$,
and the decay rates of the steepest decent contours given by integrating
$\Delta\widetilde{G}_{q}$ along $k=k_{q}-\mathrm{i}\xi$, we are required to
understand the behaviour of $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ at
$r=1-h,r,r_{0},$ and $1$ as $k\to\frac{\omega}{M},k_{r},k_{0}$, where
$\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ are given in appendix A.3.

Considering the evaluations at $r,r_{0}>1-h$ and $1$, it can be noted that we
must examine the cases that $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ are
described as $\widetilde{p}_{c1}$ and $\widetilde{p}_{c2}$ respectively (as
described in appendix A.1), or, if we have been required to perform matching,
that both are expressed in terms of $\widetilde{p}_{11}$ and
$\widetilde{p}_{12}$ (given in appendix A.2). If we examine the limit
$k\to\frac{\omega}{M}$, it will follow that in each case we are required to
take the matched solutions $\widetilde{p}_{11}$ and $\widetilde{p}_{12}$.

###  C.1 Asymptotic behaviour as $k\to\frac{\omega}{M}$

Consider first $\widetilde{p}_{c1}$ and $\widetilde{p}_{c2}$ for $k$ close to
$\frac{\omega}{M}$ and $r$ close to $1\\!-\\!h$. Since
$Q=h\sqrt{1-\omega/(Mk)}=O\big{(}(k-\frac{\omega}{M})^{\frac{1}{2}}\big{)}$,
we consider the limit $|Q|\to 0$ and set $r=1-h+RQ$ for $|R|\leq O(1)$. By
considering the recurrence formulae for the Frobenius expansion coefficients
$a_{n}$ and $b_{n}$ given in equations (56c) and (56d) in this limit, and
after some algebra, it can be found to leading order that

$\displaystyle\widetilde{p}_{c1}(1\\!-\\!h+RQ)$
$\displaystyle=Q^{3}(R-1)^{3}\left(1+\frac{3}{4}(R-1)+\frac{3}{20}(R-1)^{2}\right)+O(Q^{4}),$
(71a) $\displaystyle\widetilde{p}_{c2}(1\\!-\\!h+RQ)$
$\displaystyle=1+O(Q^{2}\log(Q)),$ (71b)
$\displaystyle\widetilde{p}_{c1}^{\prime}(1\\!-\\!h+RQ)$
$\displaystyle=3Q^{2}(R-1)^{2}\left(1+(R-1)+\frac{1}{4}(R-1)^{2}\right)+O(Q^{3}),$
(71c) $\displaystyle\widetilde{p}_{c2}^{\prime}(1\\!-\\!h+RQ)$
$\displaystyle=-Q\log(Q)\\!\left(\frac{\omega^{2}}{M^{2}}+\frac{m^{2}}{(1\\!-\\!h)^{2}}\right)\\!(R\\!-\\!1)^{2}\\!\left(\\!1+(R\\!-\\!1)+\frac{(R\\!-\\!1)^{2}}{4}\right)\\!+O(Q)$
(71d)

We consider next $\widetilde{p}_{11}$ and $\widetilde{p}_{12}$ in the same
limit. By considering the recurrence formulae for the series coefficients
$\alpha_{n}$ and $\beta_{n}$ given in equations (61a) and (61b), it can be
found that there are coefficients
$\widetilde{p}_{11}^{(n)},\widetilde{p}_{11}^{\prime(n)},\widetilde{p}_{12}^{(n)}$,
and $\widetilde{p}_{12}^{\prime(n)}$ which are $O(1)$ as $|Q|\to 0$ such that

$\displaystyle\widetilde{p}_{11}(1-h+RQ)$
$\displaystyle=\sum_{n=0}^{\infty}(RQ)^{n}\widetilde{p}_{11}^{(n)},$
$\displaystyle\widetilde{p}_{12}(1-h+RQ)$
$\displaystyle=\sum_{n=0}^{\infty}(RQ)^{n}\widetilde{p}_{12}^{(n)},$ (72a)
$\displaystyle\widetilde{p}_{11}^{\prime}(1-h+RQ)$
$\displaystyle=\sum_{n=0}^{\infty}(RQ)^{n}\widetilde{p}_{11}^{\prime(n)},$
$\displaystyle\widetilde{p}_{12}^{\prime}(1-h+RQ)$
$\displaystyle=\sum_{n=0}^{\infty}(RQ)^{n}\widetilde{p}_{12}^{\prime(n)}.$
(72b)

Note in particular that, as $|Q|\to 0$, the coefficients of the $R^{5}$ term
in both $\widetilde{p}_{11}$ and $\widetilde{p}_{12}$ tend to zero at least as
fast as $Q^{5}$, where as in $\widetilde{p}_{c1}$ from (71a) the coefficient
of $R^{5}$ tends to zero as $Q^{3}$. Hence, if we were to write
$\widetilde{p}_{c1}=A_{1}\widetilde{p}_{11}+B_{1}\widetilde{p}_{12}$, then at
least one of the coefficients $A_{1}$ and $B_{1}$ would need to tend to
infinity as $Q^{-2}$ or faster as $|Q|\to 0$. We argue that the choice of
$\widetilde{p}_{11}$ and $\widetilde{p}_{12}$ as two linearly independent
solutions about $r=1$ is arbitrary, and so by symmetry between
$\widetilde{p}_{11}$ and $\widetilde{p}_{12}$ we expect $A_{1}$ and $B_{1}$ to
be the same order of magnitude in $Q$, therefore forcing that
$A_{1}=O(Q^{-2})$ and $B_{1}=O(Q^{-2})$. Similarly, since
$\widetilde{p}_{c2}^{\prime}$ has a coefficient of $R^{4}$ which scales as
$Q\log Q$, if we were to write
$\widetilde{p}_{c2}^{\prime}=A_{2}\widetilde{p}_{11}^{\prime}+B_{2}\widetilde{p}_{12}^{\prime}$,
then at least one of the coefficients $A_{2}$ and $B_{2}$ would need to tend
to infinity as $Q^{-3}\log Q$ or faster as $|Q|\to 0$, and so we argue that
$A_{2}=O\big{(}Q^{-3}\log Q\big{)}$ and $B_{2}=O\big{(}Q^{-3}\log Q\big{)}$.

Note, however, that by evaluation the Wronskian
$\mathcal{W}(\widetilde{p}_{c1},\widetilde{p}_{c2};r)=W(r)$ at $r=1$ using
(67), and by definition of $\widetilde{p}_{11}$ and $\widetilde{p}_{12}$ at
$r=1$, considering $\mathcal{W}(\widetilde{p}_{11},\widetilde{p}_{12};r)$ at
$r=1$ shows that

$A_{1}B_{2}-A_{2}B_{1}=\frac{3(h-Q)^{2}(h+Q)^{2}(1\\!-\\!h+Q)}{4Q^{2}}=O(Q^{-2}).$
(73)

This is smaller than might have been expected from the individual scalings of
$A_{1}$, $B_{1}$, $A_{2}$ and $B_{2}$ given above, but this is expected as,
when the critical point $r_{c}^{+}$ is approached, the two linearly
independent solutions lose their linear independence, and so there is
significant cancellation between $A_{1}B_{2}$ and $A_{2}B_{1}$.

Note also from (67) that, as $|Q|\to 0$, we have

$W(r^{*})=\begin{cases}-\dfrac{3Q^{2}}{4}\\!\left(\\!1+\dfrac{Q}{1\\!-\\!h}\right)\\!&r_{0}\leq
1-h\\\ \\\
-\dfrac{3(1\\!-\\!h-r_{0})^{4}}{4r_{0}}\\!\left(\dfrac{(1\\!-\\!h)}{Q^{2}}+\dfrac{1}{Q}+O(1)\\!\right)\\!&r_{0}>1-h.\end{cases}$
(74)

Assuming that $\widetilde{p}_{11}$ and $\widetilde{p}_{12}$ are $O(1)$ when
$r$ is not close to $1-h$, it follows that

(75a) (75b) (75c) (75d) (75e) (75f)

We can use the above to establish that $\psi_{1}$ and $\psi_{2}$ are both
order 1 quantities for particular values of $r$:

$\displaystyle\psi_{1}(r)$ $\displaystyle=J_{m}(\alpha r)=O(1)$ for
$\displaystyle r<1-h;$ (75bx) $\displaystyle\psi_{2}(r)$
$\displaystyle=\widehat{C}_{2}\widetilde{p}_{1}+\widehat{D}_{2}\widetilde{p}_{2}=\widetilde{p}_{12}-\frac{\mathrm{i}\omega}{Z}\widetilde{p}_{11}=O(1)$
for $\displaystyle r>1-h.$ (75by)

We also note that $\omega-U(r^{*})k=-M(k-\frac{\omega}{M})=-\omega
Q^{2}\\!/h^{2}+O(Q^{4})=O(Q^{2})$ for $r_{0}\leq 1-h$ and is $O(1)$ for
$r_{0}>1-h$, and that

$A=-\frac{1}{3}\left(\frac{\omega^{2}}{M^{2}}+\frac{m^{2}}{r_{c}^{+2}}\right)\left(\frac{1}{Q}-\frac{1}{r_{c}^{+}}\right)-\frac{2m^{2}}{3r_{c}^{+3}}=-\frac{1}{3Q}\left(\frac{\omega^{2}}{M^{2}}+\frac{m^{2}}{r_{c}^{+2}}\right)+O(1).$
(75bz)

####  C.1.1 Behaviour of $\widetilde{G}$ as $k\to\frac{\omega}{M}$

We now use the above scalings to consider the branch point of $\widetilde{G}$
at $k=\frac{\omega}{M}$, with the aim of showing that $\widetilde{G}$ does not
experience a pole at $k=\frac{\omega}{M}$ for any value of $r_{0}$. Recall
from (32) that

$\widetilde{G}=\frac{\left(\omega-U(r^{*})k\right)}{2\upi\mathrm{i}r^{*}W(r^{*})}\frac{\psi_{1}(\check{r})\psi_{2}(\hat{r})}{C_{1}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}.$
(32)

Using the results above, if
$k=\frac{\omega}{M}+\varepsilon\mathrm{e}^{\mathrm{i}\theta}$ then for
$r_{0}\leq 1-h$,

$\widetilde{G}\sim\frac{M\varepsilon\mathrm{e}^{\mathrm{i}\theta}}{2\upi\mathrm{i}(1\\!-\\!h)}\frac{\psi_{1}(\check{r})\psi_{2}(\hat{r})}{\alpha
J_{m}^{\prime}(\alpha(1\\!-\\!h))\widehat{D}_{2}}=O(\varepsilon).$ (75ca)

If instead $r_{0}>1-h$, we find that

$\widetilde{G}\sim\frac{-Mh^{4}\varepsilon^{2}\mathrm{e}^{2\mathrm{i}\theta}\big{(}M-U(r_{0})\big{)}}{2\upi\mathrm{i}\omega(1\\!-\\!h)(1\\!-\\!h-r_{0})^{4}}\frac{\psi_{1}(\check{r})\psi_{2}(\hat{r})}{\alpha
J_{m}^{\prime}(\alpha(1-h))\widehat{D}_{2}}=O\big{(}\varepsilon^{2}\big{)}.$
(75cb)

In particular, in either case there is no pole of $\widetilde{G}$ at
$k=\frac{\omega}{M}$. Hence, we have that

$I_{\varepsilon}(x)=\frac{-1}{2\upi}\\!\\!\int_{0}^{2\upi}\\!\widetilde{G}\\!\left(\frac{\omega}{M}+\varepsilon\mathrm{e}^{\mathrm{i}\theta}\right)\exp\\!\left\\{-\mathrm{i}x\left(\frac{\omega}{M}+\varepsilon\mathrm{e}^{\mathrm{i}\theta}\right)\\!\right\\}\mathrm{i}\varepsilon\mathrm{e}^{\mathrm{i}\theta}\,\mathrm{d}\theta\to
0\qquad\text{as}\quad\varepsilon\to 0.$ (75cc)

####  C.1.2 Behaviour of $\Delta\widetilde{G}_{\frac{\omega}{M}}$ as
$k\to\frac{\omega}{M}$

We now substitute all of the above into the equation for
$\Delta\widetilde{G}_{\frac{\omega}{M}}$ given in equation (45a). First of
all, we rewrite (45a) exactly as

$\displaystyle\Delta\widetilde{G}_{\frac{\omega}{M}}=\frac{4A\big{(}\widehat{D}_{2}\psi_{1}^{-}(1\\!-\\!h)\big{)}^{2}f(r)f(r_{0})j(r^{*})}{3(1-h)Q^{3}\left(C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}+2\mathrm{i}\upi
AD_{1}\widehat{D}_{2}\right)\\!\left(C_{1}^{-}\widehat{D}_{2}-\widehat{C}_{2}D_{1}\right)},$
(75cda) $\displaystyle\text{where}\quad
j(r_{0})=-\frac{3}{4}(1-h)Q^{3}\frac{\omega-U(r_{0})k}{r_{0}W(r_{0})}\qquad\text{and}\quad
f(r)=\begin{cases}\dfrac{\psi_{1}^{-}(r)}{\psi_{1}^{-}(1\\!-\\!h)}&r<r_{0}\\\
\\\
\dfrac{D_{1}\psi_{2}^{-}(r)}{\widehat{D}_{2}\psi_{1}^{-}(1\\!-\\!h)}&r>r_{0}.\end{cases}$
(75cdb)

Taking now the leading order terms as $k\to\frac{\omega}{M}$, we find that

$\displaystyle\Delta\widetilde{G}_{\frac{\omega}{M}}$
$\displaystyle\sim-\frac{\left(\frac{\omega^{2}}{M^{2}}+\frac{m^{2}}{(1\\!-\\!h)^{2}}\right)}{4(1-h)}\frac{J_{m}\big{(}\alpha(1\\!-\\!h)\big{)}^{2}}{\alpha^{2}J_{m}^{\prime}\big{(}\alpha(1\\!-\\!h)\big{)}^{2}}f(r)f(r_{0})j(r_{0}),$
(75cea) $\displaystyle\text{where }f(r)$
$\displaystyle=\begin{cases}\dfrac{J_{m}(\alpha
r)}{J_{m}\big{(}\alpha(1\\!-\\!h)\big{)}}&r<r_{0}\\\ \\\
\dfrac{\psi_{2}^{-}(r)}{\widehat{D}_{2}}&r>r_{0}.\end{cases}$ (75ceb)
$\displaystyle\text{and }j(r_{0})$
$\displaystyle=\begin{cases}-\dfrac{\omega}{h^{2}}Q^{3}&r_{0}<1-h\\\ \\\
\dfrac{\omega(1-U(r_{0})/M)}{(r_{0}-1+h)^{4}}Q^{5}&r_{0}>1-h.\end{cases}$
(75cec)

Finally, setting $k=\frac{\omega}{M}-\mathrm{i}\xi$, so that
$Q=(1-\mathrm{i})h\sqrt{M\xi/2\omega}+O\big{(}\xi^{3/2}\big{)}$ (recalling
that $\mathrm{Re}(Q)\geq 0$), we find that $j(r_{0})$ may be written to
leading order as

$j(r_{0})=\begin{cases}\dfrac{1+\mathrm{i}}{\sqrt{2}}\dfrac{hM^{3/2}}{\omega^{1/2}}\xi^{3/2}&r_{0}<1-h\\\
\\\
-\dfrac{1-\mathrm{i}}{\sqrt{2}}\dfrac{h^{5}M^{5/2}}{\omega^{3/2}}\dfrac{1-U(r_{0})/M}{(r_{0}-1+h)^{4}}\xi^{5/2}&r_{0}>1-h.\end{cases}$
(75cf)

###  C.2 Asymptotic behaviour as $k\to k_{0}$

We now consider $k\to k_{0}$ with $r_{0}>1-h$. We have that

$\displaystyle r_{0}-r_{c}^{+}$ $\displaystyle=-\frac{\omega
h^{2}(k-k_{0})}{2Mk_{0}^{2}Q_{0}}+O\big{(}(k-k_{0})^{2}\big{)},$
$\displaystyle\text{where}\qquad Q_{0}=h\sqrt{1-\frac{\omega}{Mk_{0}}}.$
(75cg)

Hence, in this limit, $\widetilde{p}_{1}(r_{0})$ and
$\widetilde{p}_{2}(r_{0})$ may always be evaluated in terms of
$\widetilde{p}_{c1}$ and $\widetilde{p}_{c2}$, as we are always eventually
within their radius of convergence. Hence, in this limit,

$\displaystyle\widetilde{p}_{1}(r_{0})$ $\displaystyle=\left(\frac{-\omega
h^{2}(k-k_{0})}{2Mk_{0}^{2}Q_{0}}\right)^{\\!3}+O\big{(}(k-k_{0})^{4}\big{)},$
$\displaystyle\widetilde{p}_{2}(r_{0})$
$\displaystyle=1+O\big{(}(k-k_{0})^{2}\big{)}.$ (75ch)

For $r\not=r_{0}$, the Bessel function, Hankel functions, and
$\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ all behave as $O(1)$ quantities
when evaluated at $1-h,r$, and $1$, resulting in $O(1)$ behaviour for
$C_{1},D_{1},\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{C}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{C}{\scalebox{-1.0}{\tmpbox}}_{2},\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{D}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{D}{\scalebox{-1.0}{\tmpbox}}_{2},\widehat{C}_{2}$
and $\widehat{D}_{2}$. It can be shown that $A=O(1)$ and that

$W(r_{0})=-\frac{3h^{4}\omega^{2}(k-k_{0})^{2}}{4Q_{0}^{2}M^{2}k_{0}^{4}}+O\big{(}(k-k_{0})^{3}\big{)}.$
(75ci)

####  C.2.1 Behaviour of $\widetilde{G}$ as $k\to k_{0}$ and the residue of
the non-modal $k_{0}$ pole

Substituting all the above into (32) (as $k\to k_{0}$ from above) gives

$\widetilde{G}(k)=\frac{-2Mk_{0}^{2}(\omega-
Mk_{0})}{3\upi\mathrm{i}r_{0}h^{2}\omega(k-k_{0})}\frac{1}{C_{1}^{+}\widehat{D}_{2}-\widehat{C}_{2}D_{1}}\begin{cases}\widehat{D}_{2}\psi_{1}(r)&r<r_{0}\\\
D_{1}\psi_{2}(r)&r>r_{0}\end{cases}+O(1),$ (75cj)

confirming a pole at $k=k_{0}$ that gives a residue contribution once
integrated around of

$R_{0}^{+}(k_{0})=\frac{2Mk_{0}^{2}(\omega-
Mk_{0})\mathrm{e}^{-\mathrm{i}k_{0}x}}{3\pi
r_{0}h^{2}\omega(C_{1}^{+}\widehat{D}_{2}-\widehat{C}_{2}D_{1})}\begin{cases}\widehat{D}_{2}\psi_{1}(r)&r<r_{0}\\\
D_{1}\psi_{2}(r)&r>r_{0}.\end{cases}$ (75ck)

####  C.2.2 Behaviour of $\Delta\widetilde{G}_{0}$ as $k\to k_{0}$

Moreover, we may substitute all the above into $\Delta\widetilde{G}_{0}$ from