The critical layer in quadratic flow boundary layers over acoustic linings

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 Green's function is constructed using Frobenius series. The critical layer gives a non-modal contribution from the continuous spectrum branch cut, and dominates the downstream pressure perturbation in certain cases, particularly for thicker boundary layers. The continuous spectrum branch cut is also found to stabilize what are otherwise convectively unstable modes by hiding them behind the branch cut. Overall, the contribution from the critical layer is found to give a neutrally stable non-modal wave when the source is located within the sheared flow region, and to decay algebraically along the duct as $O(x^{-{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 lose accuracy, is also able to locate modal poles hidden behind the branch cut, which other methods are unable to find; this includes the stabilized hydrodynamic instability. Matlab code is provided to compute the Green's 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 engine noise. With aircraft noise being subjected to ever increasing restrictions, being able to model this noise successfully becomes increasingly important. In particular, aircraft engine noise at take-off 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 ${\rm Re}(p\exp \{\mathrm {i}\omega t\})$ leads to an oscillating normal boundary velocity ${\rm 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; Khamis & Brambley 2017; Schulz et al. 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 consider only 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 time step owing to the Courant–Friedrichs–Lewy (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 (2.5), a second-order linear ordinary differential equation for the pressure perturbation within the duct due to Pridmore-Brown (1958). The Pridmore-Brown equation has been the subject of much analysis (e.g. Mungur & Gladwell 1969; Ko 1972; Swinbanks 1975; Nagel & Brand 1982; Brambley, Darau & Rienstra 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 used frequently, 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, McAlpine & Astley 2010; Oppeneer, Rienstra & Sijtsma 2016; Rienstra 2021).

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.

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 (2.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 of 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, Davis & Peake 2012b; Posson & Peake 2013; Mathews & Peake 2018b). For this reason, the Green's function 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 of 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, Darau & Brambley 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 that 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 non-zero second derivative of the mean flow profile, $U''(r)$, but for the constant-then-linear mean flow, $U''(r)$ either is 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 no-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''(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 § 2, along with a derivation of the Pridmore-Brown equation by spatially Fourier transforming the linearized Euler equations. The Frobenius expansion is then used in § 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 § 4 by evaluating numerically 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 § 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$:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2. Homogeneous solutions within the region of uniform flow

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

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

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

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

it will turn out later that the branch chosen for $\alpha$ does not matter, although for definiteness, one may choose $\textrm {Re}(\alpha )>0$. Bessel's equation has two pairs of linearly independent solutions that we will make use of: the Bessel functions of the first and second kind, $J_m(\alpha r)$ and $Y_m(\alpha r)$; and the Hankel functions of the first and second kind, $H_m^{(1)}(\alpha r)$ and $H_m^{(2)}(\alpha r)$. More information regarding these can be found in Abramowitz & Stegun (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 of the homogeneous Pridmore-Brown equation (2.9) when $U(r)$ varies by proposing a Frobenius expansion about the singularities of the Pridmore-Brown equation.

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

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

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

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

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

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

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

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

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

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

and that

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

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

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

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

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

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

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

2.4. Homogeneous solutions across the full domain

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

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

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

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

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

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

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

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

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

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

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

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

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

where the factor at the beginning comes from the Wronskian of $H_m^{(1)}$ and $H_m^{(2)}$ from Abramowitz & Stegun (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 $\tilde {p}_2(r)$ when $r< r_c^{+}$, we have, provided that $r_c^{+}<1$,

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

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

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

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

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

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

2.5. Modal solutions

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

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

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

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

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

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

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

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

where we have made use of the Bessel function identities 9.1.3, 9.1.4 and 9.1.16 from Abramowitz & Stegun (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$, because $\psi _1$ and $\psi _2$ are both $C^{1}$ continuous, it follows that for $0\leq r\leq 1-h$, we can set $r\,\mathcal {W}(\psi _1,\psi _2;r)=(1-h)\,\mathcal {W}(\psi _1,\psi _2;1-h)$. We therefore arrive at the conclusion that

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

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

3. Inhomogeneous solutions and inverting the Fourier transform

3.1. Inhomogeneous solution of the Pridmore-Brown equation

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

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

where

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

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

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

3.2. Analytic continuation behind the critical layer branch cut

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

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

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

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

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

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

Similarly, we may rewrite (3.5) as

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

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

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

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

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

3.3. Inverting the Fourier transform

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

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

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

3.3.1. Choosing an inversion contour

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

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

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

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

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

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

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

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

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

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

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

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

3.3.3. Contribution from the critical layer branch cut

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

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

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

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

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

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

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

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

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

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

The steepest descent integrals are defined as

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

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

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

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

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

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

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

The critical layer branch cut contribution (3.13) contains integrals $I_q(x)$ given by (3.15) that are amenable to asymptotic analysis in the limit $x \to \infty$, using the method of steepest descent. Having already deformed the integration contours onto the steepest descent contours, so that the integrands have had their $x$-dependent oscillation removed and are now exponentially decaying along the contour, we may directly apply Watson's lemma (Watson 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$,

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

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

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

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

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

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

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

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

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

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

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

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

3.4.1. Comparisons with previous far-field scalings

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

Table 1. Comparison of the different decay rates given for a general flow profile by Swinbanks (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 difference in the behaviour of the $I_{{{\omega }/{M}}}$ integrals may be understood as having two causes. The first is the difference in behaviour of the constant $A$, given in general as

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

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

4. Numerical results

In this section, the above analysis is illustrated with some numerical examples. The Frobenius series solutions $\tilde {p}_1$ and $\tilde {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, using the Matlab code in the supplementary material. In particular, this is more expensive numerically and prone to numerical rounding errors near the edge of the radius of convergence for each of the solutions, and care must therefore be taken to choose an accurate and efficient numerical implementation of $\tilde {p}_1$ and $\tilde {p}_2$ in terms of the Frobenius series solutions. 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 $|\textrm {Im} (k)|<400$ have been included.

Throughout this section, we show results from four parameter sets, detailed in table 2. 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 take-off, while parameter set C might represent the same type of mode during the landing approach.

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

In § 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 § 4.2, we compare the various contributions from the critical layer branch cut described in § 3.3.3, including their large-$x$ behaviour, and show that these agree with the predicted large-$x$ behaviour from § 3.4. In § 4.3, the full solution in terms of $x$ and $r$ is plotted, and these results are compared to the linear boundary layer case. Finally, in § 4.4, we investigate how the results vary as we vary parameters, including the frequency $\omega$, the boundary layer thickness $h$, the wall impedance $Z$, and the steady mean flow velocity $M$.

4.1. Pole locations

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

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

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

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

4.2. Branch cut contributions

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

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

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

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

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

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

4.3. Full Fourier inversion

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

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

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

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

Figure 11 compares the wave field generated in a quadratic boundary layer with the wave field in a linear boundary layer profile (as studied by Brambley et al. 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 found to be unstable only if the boundary layer is sufficiently thin. Even when both flow profiles give rise to a convective instability, we can see in figure 11 that the growth rates of the instabilities can be significantly different.

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

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

4.4. Variation of results with changing parameters

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

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

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

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 no-slip boundary condition at the duct wall, whilst maintaining a $C^{1}$ continuous flow. For such a flow profile, irrespective of the the thickness of the boundary layer, a solution of the Pridmore-Brown equation has been constructed making use of two Frobenius series expansions, valid for any wavenumber $k$. This enables the evaluation of the Green's function of the Pridmore-Brown equation, which 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.

In this work, we have aimed to construct the Green's function solutions, which is equivalent to the solution for a point mass source in the linearized Euler equations. The Green's function is in some sense the general solution, as the solution subject to any forcing can be written as an integral over the Green's function, suitably weighted. Because of this, any behaviour the equations are capable of must necessarily be demonstrated in the Green's function solution, so once the behaviour of the Green's function is understood, the equations cannot hold any further surprises. This is particularly important in this case, considering that the Green's function solution has been shown to include non-modal contributions, such as the critical layer branch cut, that cannot be investigated clearly using other methods, such as eigenfunction methods that capture only the modal contributions.

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, whereas 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 that 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 with 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 summarized in § 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, present only when $r_0>1-h$, which does not decay with distance from the point source and is therefore dominant in the far field 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, travelling 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_{{{\omega }/{M}}}$, $I_r$ and $I_0$, which are the result of accounting for the branch points coming from the critical points of the Pridmore-Brown equation. These contributions have phase speeds equal to, respectively, the uniform flow speed $M$, $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 flow speed is sufficiently slow, although the specific stability boundary also depends on the impedance $Z$ and frequency $\omega$.

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

The various decay rates of the components of the critical layer have been predicted previously by Swinbanks (1975) and Brambley et al. (2012a); a summary can be found in table 1. Swinbanks (1975) considered only the contribution from waves with phase velocity $U(r_0)$, which are present only 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 (1975) considering a two-dimensional flow in a rectilinear duct with an arbitrary flow profile, and Brambley et al. (2012a) considering only a constant-then-linear flow profile in a three-dimensional cylindrical duct; in particular, Swinbanks (1975) emphasized 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 as 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^{-{5}/{2}})$ is found, while when the source is within the sheared flow, we instead have a faster rate of decay of $O(x^{-{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 that 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 a 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 predicted accurately 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$, so $k_0$ and $k^{+}$ coinciding could be expected to occur for any parameters leading to exact neutral stability. One could also extend this problem to a non-constant mean density and sound speed making use of (2.5). For a given mean density profile $\rho _0(r)$ and sound speed profile $c_0(r)$, one could still construct a solution of the resulting Pridmore-Brown equation, taking careful notice of the potentially complex roots of $c_0(r)$. It would still be possible to construct a solution using the Frobenius series solutions so long as these are not double roots or of higher order, and ${\rho _0^{\prime }(r)}/{\rho (r)}$ has at most regular singularities in the complex $r$-plane. The critical layer branch cut would still occur in identical form to that seen in our work (except in the case where these singularities occur at $r=r_c^+$), although the resulting scaling in the various limits seen above may vary. When calculating the scalings in the various limits above, the work given here would provide a suitable outline for the approach to be taken. 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).