The critical layer in quadratic flow boundary layers over acoustic linings

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


Introduction
The propagation of sound through an otherwise steady mean flow has many important applications. One such application is predicting and optimizing aircraft engines noise. With aircraft noise being subjected to ever increasing restrictions, being able to successfully model this noise becomes increasingly important. In particular, aircraft engine noise at takeoff depends critically on the sound absorbing performance of acoustic liners. Unfortunately, acoustic liner performance in the presence of a steady mean flow is poorly predicted by existing theory, as demonstrated by comparisons to laboratory experiments (Renou & Aurégan 2011;Spillere, Bonomo, Cordioli & Brambley 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 Re(p exp{iωt}) leads to an oscillating normal boundary velocity Re(v exp{iωt}) given by p = Z(ω)v. Such impedance boundary conditions are well understood for a mean flow that satisfies no-slip at the boundary. Often, however, we use a simplified model where the mean flow does not satisfy no-slip at the boundary: for example, uniform axial flow in a duct. For slipping mean flows, it is known that the impedance boundary condition must be modified. A common modified boundary condition is the Myers, or Ingard-Myers, boundary condition (Ingard 1959;Myers 1980). This boundary condition is known to be the correct limiting behaviour for an inviscid mean flow boundary layer in the limit that the boundary layer thickness tends to zero (Eversman & Beckemeyer 1972;Tester 1973). However, this boundary condition, when applied in the time domain, is ill-posed (Brambley 2009). Several alternative boundary conditions have been suggested (Brambley 2011b;Schulz, Weng, Bake, Enghardt & Ronneberger 2017;Khamis & Brambley 2017;Aurégan 2018), which each attempt to include more relevant physics, including the effect of the mean flow boundary layer and the effect of viscosity. However, these boundary conditions come with their own complications, including the need to fit further free parameters, and as yet none have been made to agree with laboratory experiments (Spillere et al. 2020).
In light of this difficulty with boundary conditions in slipping mean flow, one may instead only consider mean flows U (r) that satisfy no-slip at the boundary (e.g. Weng, Schulz, Ronneberger, Enghardt & Bake 2017). Doing so, however, involves solving for the sound in a strongly varying mean flow, which is especially taxing when the boundary layers are particularly thin. Numerically resolving the sound in thin boundary layers requires a fine resolution, which then also requires a small timestep owing to the CFL condition. Progress may be made analytically by considering the simplified situation of a straight rectilinear or cylindrical duct containing axial mean flow (as depicted later in figure 1). By linearizing the Euler Equations about this steady mean flow and assuming exp{iωt−ikx} dependence, one eventually arrives at the Pridmore-Brown equation (2.5), a secondorder linear ODE for the pressure perturbation within the duct due to Pridmore-Brown (1958). The Pridmore-Brown equation has been the subjected of much analysis (e.g. Mungur & Gladwell 1969;Ko 1972;Swinbanks 1975;Nagel & Brand 1982;Brambley, Darau & Rienstra 2012a;Rienstra 2020), owing to its complexity. One complexity is that, treating the frequency ω as known and solving for the axial wavenumber k as the eigenvalue, the Pridmore-Brown equation is not Sturm-Liouville and results in a nonlinear eigenvalue problem for k. A second complexity is that the Pridmore-Brown equation possesses a regular singularity, referred to as a critical layer or continuous spectrum. Despite these difficulties, eigenfunction expansions using eigenfunctions of the Pridmore-Brown equation are frequently used, with the eigenfunctions assumed to form a complete basis (despite the problem being non-self-adjoint) and the effect of the critical layer ignored (e.g. Brooks & McAlpine 2007;Olivieri, McAlpine & Astley 2010;Oppeneer, Rienstra & Sijtsma 2016;Rienstra 2021).
The lack of completeness of the modal solutions of the Pridmore-Brown equation motivates the investigation of the Green's function solution. The Green's function is the solution of the governing equations subject to a point forcing; for example, a point mass source leads to the right-hand-side of equation (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 to the governing equations, in contrast to a modal eigenvalue expansion which can only express an arbitrary solution if the modal basis is complete. The Green's function is also worth considering on its own merits without reference to a particular forcing, since if the governing equations are capable of exhibiting a particular feature (such as instability, focusing, perfect reflection, etc), then the Green's function must also exhibit that feature. The Green's function is also used in various approximation techniques (e.g. Brambley, Davis & Peake 2012b;Posson & Peake 2013;Mathews & Peake 2018b). For this reason, the Green's functions has been constructed for a variety of acoustical situations (e.g. Rienstra & Tester 2008;Brambley et al. 2012a;Mathews & Peake 2017, 2018a. In particular, the Green's function solution to the Pridmore-Brown equation naturally includes the critical layer. The critical layer, or continuous spectrum, is a singularity of the linearized Euler equations occurring when the phase velocity of the perturbation, ω/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 the pressure perturbation from the critical layer at a radius r consisted of three distinct components with phase velocities U (0), U (r) and U (r 0 ), each with different decay rates. However, Brambley et al. (2012a) chose a rather special mean flow profile. In particular, the critical layer is usually caused by a nonzero second derivative of the mean flow profile, U (r), but for the constant-then-linear mean flow U (r) is either identically zero or has a delta function discontinuity; in the constant-then-linear case, Brambley et al. (2012a) instead attributed the critical layer to the cylindrical geometry.
In many cases, the effect of the critical layer is negligible in comparison with the modal sum of the acoustics modes. However, when all acoustic modes are cut-off and nonpropagating, the effect of the critical layer will be dominant. Moreover, Brambley (2013, figure 6) showed that a mode representing a hydrodynamic instability could interact with the critical layer, although this was not seen for a constant-then-linear mean flow profile.
Since the critical layer is a singularity of the Pridmore-Brown equation, traditional numerical methods are particularly inaccurate near the critical layer. This often manifests as a collection of spurious numerical modes being located along the critical layer. In contrast, previous studies have used a Frobenius expansion about the singular point r = r c (e.g. Heaton & Peake 2006;Campos & Kobayashi 2009;Brambley et al. 2012a). This technique both gives increasing accuracy as the critical layer is approached, and allows analytical continuation behind the critical layer branch cut. For example, Brambley et al. (2012a, figure 10) found a previously unknown mode close to the critical layer that was unable to be resolved numerically using more traditional finite differences. One complication of the Frobenius series, however, is that, much like a power series, it has an associated radius of convergence. For the constant-then-linear mean flow Frobenius expansion (Brambley et al. 2012a), this did not prove a problem, as the radius of convergence covered the region of interest in all cases that were considered. For general flow profiles this will not be the case, and a solution covering the entire region of interest will involve multiple Frobenius expansions with overlapping radii of convergence; this will turn out to be the case here. By matching two different expansions in a region where both converge, a hybrid solution may be constructed that spans the whole region of interest.
Here, we use the Frobenius expansion method as described by Brambley et al. (2012a), and apply it to a mean flow that is constant in the centre of the duct and then varies quadratically within a boundary layer to satisfy non-slip at the wall. As well as being more realistic than the constant-then-linear profile considered by Brambley et al. (2012a), this mean flow profile is twice differentiable, allowing U (r) to enter the analysis, and as such we expect the results to be more representative of an arbitrary mean flow profile. The Frobenius expansion is derived in section 2, along with a derivation of the Pridmore-Brown equations by spatially Fourier transforming the linearized Euler equations. The Frobenius expansion is then used in section 3 to derive the Green's function for a point mass source, including inverting the spatial Fourier transform and investigating the farfield behaviour. Results are presented in section 4 by numerically evaluating the Frobenius expansions and the Green's function. These results are compared against previous results, particularly against the predictions by Swinbanks (1975) and the constant-then-linear results by Brambley et al. (2012a). Finally, the implications of this work are discussed, and areas for further research highlighted, in section 5.

Constructing the Pridmore-Brown Equation
The governing equations for what follows are the Euler equations with a mass source q, where u is the fluid velocity, p is the pressure, ρ is the density, and c 2 = ∂p ∂ρ | s is the square of the sound speed. In what follows, we take the mass source q to be a small timeharmonic point mass source. In cylindrical coordinates (x, r, θ), with a suitable choice of origin, this mass source q may in general be taken as where is the small amplitude, ω 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 , 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 only needed to leading order in what follows. Without loss of generality, all perturbations are expanded using a Fourier series in θ and a Fourier Transform in x. As a result, the pressure perturbation is given aŝ (2.4) and similarly for the densityρ and the velocity componentsû,v andŵ. Substituting these into the Euler equations (2.1), and linearizing by ignoring terms of O( 2 ) or smaller, each ofρ,ũ,w, and finallyṽ may be eliminated, to leave a second order ODE in the radial coordinate r for p, where a prime denotes the derivative with respect to r. This is the Pridmore-Brown (1958) equation for a point mass source, written in cylindrical co-ordinates.
One boundary condition to (2.5) is regularity at r = 0. The singular solution behaves, for m = 0, as O(r −|m| ) as r → 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 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, where Z(ω) is the impedance of the duct wall, and the two expressions are equivalent in light of (2.5b). A hard wall corresponds to Z → ∞, and hence toṽ(a) = 0, or equivalently top (a) = 0.
In what follows, we make the simplifying assumption of a constant density ρ 0 (r). This Figure 1. A cross sectional view of a cylindrical duct with lined walls containing sheared axial flow. ρ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.
is a homentropic assumption, and implies that c 0 (r) is also constant. We may then nondimensionalize speeds by the sound speed c 0 , densities by ρ 0 , and distances by the duct radius a. Note that this places the impedance boundary condition in nondimensional terms at r = 1. We also assume a flow profile U (r) that is uniform, except within a boundary layer of width h where it varies quadratically: With the nondimensionalization 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 to the homogeneous form (2.9)

Homogeneous Solutions Within the Region of Uniform Flow
Within the region of uniform flow, the homogeneous Pridmore-Brown equation (2.9) reduces to This is Bessel's equations of order m rescaled by α, where (2.11) it will turn out later that the branch chosen for α does not matter, although for definiteness one may choose Re(α) > 0. Bessel's equation has two pairs of linearly independent solutions that we shall make use of: the Bessel functions of the First and Second kind, J m (αr) and Y m (αr); and the Hankel functions of the first and second kind. H (1) m (αr) and H (2) m (αr). More information regarding these can be found in Abramowitz & Stegun (1964). It is worth noting that only J m (αr) is regular at r = 0, with the other solutions all requiring a branch cut along αr < 0, with a singularity at αr = 0.

Homogeneous Solutions Within the Region of Sheared Flow
In this section, we will construct the solution to 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 ω − 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 ω − 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) For convenience, we will take Re(Q) 0, so that Re(r + c ) 1 − h and Re(r − c ) 1 − h. Since solutions with this quadratic flow profile U (r) are only valid for 1 − h < r < 1, it will therefore be r + c that we are mostly concerned about here. Following Brambley et al. (2012a), we propose a Frobenius expansion (Teschl 2012) about the regular singularity r + c , Specifying that a 0 = 0 results in a condition on σ, and we find that σ = 0, 3. By Fuchs theorem (Teschl 2012), this gives a pair of linearly independent solutions of the form The coefficients a n and b n are derived in appendix A, where, in particular, it is found that and that (2.16) the latter in agreement with equations (2.3)-(2.5) of Brambley et al. (2012a). We note in passing that in practice we may be limited by the radius of convergence of (2.14), and in such cases the solutions given above are analytically continued by a companion expansion of the Pridmore-Brown equations about r = 1, as described in appendix A.2.
Other than being a complication concerning numerical convergence, this complication may be ignored, and p 1 and p 2 thought of as being defined by the expressions in (2.14).
Due to the log term in 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 ∈ [1 − h, 1], and so the branch cut must avoid crossing the real r axis between 1 − h and 1. In the following, we achieve this by choosing the branch cuts parallel to the imaginary axis and away from the real axis, as depicted in figure 2. When r + c is real and 1 − h < r + c < 1, no suitable choice of branch cut exists, and as a result any solution p(r) with p(r + c ) = 0 necessarily has a singular third derivative at r + c . This only occurs for particular values of k, however, and we can map the corresponding values of k in the complex k plane to find they fall exactly on the half line [ ω M , ∞); this range of excluded values of k we refer to as the critical layer branch cut. As r + c becomes real, note that the value of p 2 (r + c ) is different depending on whether we approach from positive or negative imaginary part. Thinking of r + c (k) as a function of k, this corresponds to approaching the critical layer branch cut [ ω M , ∞) in k from above or below. This re-enforces the consideration of the critical layer appearing as a branch cut in the complex k plane. The change in p 2 when crossing the critical layer branch cut from below to above is described as Where H(r) is the Heaviside function.
In order to retrieve this result we need only consider the log(r − r + c ) term of p 2 . Note that ∂r + c /∂k > 0 for real k and real positive ω; hence, if k is nearly real and Im(k) > 0, then Im(r + c ) > 0, and we must take the branch cut of log(r − r + c ) upwards towards +i∞. Similarly, for Im(k) < 0 then Im(r + c ) < 0, and the branch cut for log(r − r + c ) must be taken downwards to −i∞. 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 −π. When we choose the branch cut into the lower half plane, this corresponds to a complex argument of π. This difference results in the jump of 2πi given. If we instead consider r > r + c , the same argument is retrieved regardless of which direction we take the branch-cuts, and so no jump is observed. This is the reason for the presence of the Heaviside function.

Homogeneous Solutions Across the Full Domain
In order to construct a full solution in r ∈ [0, 1], we now construct two solutions ψ 1 (r) and ψ 2 (r) that solve (2.9) across r ∈ [0, 1], by patching together the solutions derived above in sections 2.2 and 2.3. We construct ψ 1 (r) to satisfy the boundary condition at r = 0 (2.6), by taking where the coefficients C 1 and D 1 ensure C 1 continuity, and are given by 19b) and W (r) = W( p 1 , p 2 ; r) is the Wronskian of p 1 and p 2 , given in appendix A.4 as Having constructed ψ 1 to satisfy the boundary condition at r = 0, we now proceed to construct ψ 2 which satisfies the boundary condition (2.7) at r = 1. Writing ψ 2 in terms of the homogeneous solutions derived above, we choose C 2 and D 2 to satisfy ψ 2 (1) = 1 and ψ 2 (1) = − iω Z . This forces a non-zero normalized solution to ψ 2 which satisfies the boundary condition (2.7) at r = 1, and leads to . (2.22) The coefficients C 2 and D 2 are chosen such that our solution is C 1 continuous at r = 1−h, giving where the factor at the beginning comes from the Wronskian of H We will also require later the jump in behaviour of ψ 1 and ψ 2 as k crosses the critical layer branch cut from below to above. Since any jump comes from the log term in p 2 (r) when r < r + c , we have, provided r + c < 1, Note that, if r + c > 1, then ∆ψ 1 = ∆ψ 2 = 0, since the ψ 1 and ψ 2 solutions are uniquely defined by their boundary conditions and no branch point occurs on the interval r ∈ [1 − h, 1] to cause a jump.

Modal solutions
Modal solutions of the homogeneous Pridmore-Brown equation (2.9) are nonzero solutions p(r) that satisfy both the boundary conditions at r = 0 and at r = 1 (2.6,2.7). In general, satisfying both boundary conditions would force the solution p(r) ≡ 0, so nonzero solutions exist only for particular modal eigenvalues k (assuming ω is given and fixed). In contrast, the solution ψ 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 ψ 1 (r). Likewise, the solution ψ 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 ψ 2 (r). In general, ψ 1 and ψ 2 are linearly independent, and so their Wronskian (2.26) However, if p(r) is nonzero and satisfies both boundary conditions at r = 0 and r = 1, then p(r) = aψ 1 (r) = bψ 2 (r) for some nonzero coefficients a, b. In other words, a modal solution is one where ψ 1 and ψ 2 are linearly dependent, and so W(ψ 1 , ψ 2 ; r) ≡ 0.
For 1 − h r 1, substituting ψ 1 from (2.18) and ψ 2 from (2.21) into the Wronskian (2.26) gives where W (r) is the Wronskian between p 1 and p 2 and is given earlier in (2.20). Since p 1 and 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) = 0 except at the critical layer r = r + c . A modal solution, therefore, is given by the condition that C 1 D 2 − C 2 D 1 = 0, which is independent of r, and implies that C 1 /D 1 = C 2 / D 2 and so that ψ 1 and ψ 2 are multiples of one another.
The same can be seen for r 1 − h. In this case, the Wronskian (2.26) becomes 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 rW(ψ 1 , ψ 2 ; r) is a constant independent of r for 0 r 1 − h. Since W(ψ 1 , ψ 2 ; r) is continuous in r across r = 1 − h, since ψ 1 and ψ 2 are both C 1 continuous, it follows that for 0 r 1 − h we can set rW(ψ 1 , ψ 2 ; r) = (1 − h)W(ψ 1 , ψ 2 ; 1 − h). We therefore arrive at the conclusion that and that a mode corresponds to the dispersion relation 0 = D(k, ω) = C 1 D 2 − C 2 D 1 .
In the next section, we see how these modal solutions occur naturally as poles in the solution of the non-homogeneous Pridmore-Brown equation.

Inhomogeneous Solution to the Pridmore-Brown Equation
While previously we have only been solving the homogeneous form (2.9), our original problem was to solve the inhomogeneous Pridmore-Brown equation (2.5a) subjected 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 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, G may be written as a multiple of the homogeneous solution ψ 1 for r < r 0 and as a multiple of the homogeneous solution ψ 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 and once again W(ψ 1 , ψ 2 ; r) is the Wronskian of ψ 1 and ψ 2 . Using (2.29), this may be rewritten as

Analytic continuation behind the critical layer branch cut
The solution for G in (3.3) above contains a branch cut along the critical layer k ∈ [ ω M , ∞). We now introduce the following additional notation. When evaluating a function f (k) on the branch cut, Note that the definition of ∆f agrees with the use of ∆ in equations (2.17,2.24,2.25) above. By using these equations, we find that A typical branch cut, such as the branch cut in √ 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 ∈ [1 − h, 1]. None-the-less, noting from (3.5) that ∆ G is well defined function for general complex k, we may use equation (3.5) to analytically continue G behind the critical layer branch cut. For real ω, we therefore define the analytic continuation of G behind the branch cut into the lower-half k-plane as (3.6) Similarly, we may rewrite (3.5) as which allows the analytic continuation of G into the upper-half k-plane, (3.8) The utility of these analytic continuations in not readily apparent. However, their use allows for poles of G, corresponding to modal solutions to the homogeneous Pridmore-Brown equation, to be tracked behind the branch cut, and in particular a possible hydrodynamic instability mode will later be found to be hidden behind the critical layer branch cut in certain cases. Their use also allows the deformation of integral contours behind the critical layer branch cut, as will be needed for the steepest descent contours needed for the large-x asymptotic evaluation of the inverse Fourier transform.
In what follows k + and k − denote modal poles, see section 2.5, of only G + or G − respectively.

Inverting the Fourier Transform
Having formulated G as the solution to the inhomogeneous Pridmore-Brown equation (2.5a), to recover the actual pressure perturbationp(x, r, θ), 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 Note, however, that the critical layer branch cut is located along the real-k axis k ∈ [ ω M , ∞). We are therefore required to be careful in choosing a suitable inversion contour C.

Choosing an Inversion Contour
In order to choose the correct Fourier inversion contour 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 disturbancep), which is otherwise lost when considering a time-harmonic forcing, as we do here. A more in-depth description is available in many places in the literature (e.g. Brambley 2009, appendix A).
In order to make use to the Briggs-Bers criterion, the rate of exponential growth of the solution must be bounded; that is, there must exist Ω, K > 0 such that, if Im(ω) < −Ω, then G is analytic for |Im(k)| < K. For a given ω with Im(ω) < −Ω, we take the k-inversion contour 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 ω that are of interest, the imaginary part of ω is smoothly increased to 0, and the locations of any singularities tracked throughout this process. During this process, the k-inversion contour C must be smoothly deformed in order to maintain analyticity; that is, no singularities must cross the k-inversion contour. Assuming this process may be completed and Im(ω) increased to zero, then the resulting k-inversion contour C is the correct causal contour. Since for x < 0 the exp{−ikx} term is exponentially small as |k| → ∞ in the upper-half k-plane, for x < 0 we may close the contour with a large semi-circular arc at infinity in the upper-half k-plane, denoted C > . The resulting contours, for real ω, are illustrated in 3, for a typical unstable case. The majority of singularities of G are poles which do not cross the real k axis as Im (ω) 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 Im(ω) 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 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 for some critical radius r c , and so is found in the lower-half k-plane for Im(ω) < 0. Thus, as shown in figure 3, for x > 0 in order to close C in the lower-half k-plane, we must pass around the critical layer branch cut, denoted by the contour C b , before closing in the lower-half k-plane with a semi-circular arc denoted C > . The contribution from integrating around the critical layer branch cut, C b , leads to the non-modal contribution of the critical layer, and is discussed in detail below in section 3.3.3.

Contribution from the poles of G
We may now write the integral around the closed contour as a sum of residues of poles: where I(x) is the contribution from integrating around the critical layer branch cut contour C b discussed in the next section, R(k j ) is the residue from a pole at k j discussed below, and the notation Im(k j ) > Im(C) is used to denote poles k j lying above the inversion contour C.
The poles of G correspond to zeros of the denominator of 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 D 2 − C 2 D 1 = 0. As discussed in section 2.5, this occurs when both ψ 1 and ψ 2 satisfy both boundary conditions at r = 0 and r = 1. These modal poles can be further classified into acoustic modes and surface modes: acoustic modes are those for which α in equation (2.11) has a small imaginary part, and correspond to functions which are oscillatory in r; and surface modes are those for which α has a significant imaginary part, and correspond to functions which decay exponentially away from the duct walls at r = 1. For different parameters, we may find a variety of surface modes, and two with which we will be particularly interested here will be denoted k − and k + . For further details of surface modes, the reader is referred to the existing literature (e.g. Rienstra 2003;Brambley 2013).
Since the modal poles occur as zeros of C 1 D 2 − C 2 D 1 = 0, which we shall assume are simple zeros, the contribution from the residues of these poles are given by (3.11) The second type of poles are the non-modal poles, which occur when W (r * ) = 0. These occur when we loose independence between p 1 and p 2 at r * . Note from the formula for W (r), equation (2.20), that W (r * ) = 0 implies that r * = r + c . Since 1 − h r * 1 this can only occur when k is located on the critical layer branch cut. In what follows, we will refer to this non-modal pole as k 0 . Note that k 0 is a function of the radial location of the point source r 0 (through r * ), which is unlike the modal poles for which k j is independent of the value of r 0 ; this is one reason this k 0 pole is referred to as a nonmodal pole. However, since our closed contour goes around the critical layer branch cut (along contour C b ), this pole is always excluded from the sum of residues in (3.10) above, and only occurs within the calculation of I(x), which we consider next.

Contribution from the critical layer Branch Cut
The contribution from the critical layer branch cut, including any non-modal pole k 0 along the branch cut, is contained solely within the integral along around the critical layer branch cut denoted C b in figure 3, However, as it stands, this integral for I(x) is oscillatory, owing to the e −ikx factor in the integrand, and so is difficult to accurately compute numerically. This is especially true for large values of x. Instead, it is helpful to deform the integral onto the Steepest Descent contour, for which e −ikx is exponentially decaying along the contour. This contour deformation has three benefits: firstly, it allows accurate numerical calculation of the integral; secondly, it allows the derivation of large-x asymptotics using the Method of Steepest Descents; and thirdly, it brings insight into the various contributions that make up I(x). In deforming the integration contour, however, we must analytically continue G behind the branch cut (as described in section 3.2 above), and carefully deform around any poles and branch points. This is illustrated schematically in figure 4. Note that poles and branch points of G may exist behind the critical layer branch cut, and we must therefore use analytic continuations of G; the reader is reminded that G + is the analytic continuation of G down behind the branch cut from above, while G − is the analytic continuation of G up behind the branch cut from below. Here, we use the notation that a pole of G + with Re(k) > ω M is denoted k + , and a pole of G − with Re(k) > ω M is denoted k − . Thus, a k + pole with Im(k + ) < 0 or a k − pole with 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 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 ω M depicted in figure 4, depending on the location of the k + pole.
The Steepest Descent contours are where e −ikx is exponentially decaying; i.e. towards −i∞ in the complex k plane. There is no difficulty deforming the contour at infinity, since e −ikx is exponentially small there (provided x > 0, which is the only case in which the critical layer branch cut contributes). Along the branch cut there are up to three branch points singularities, denoted ω 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 p 2 (r), and the presence of p 2 (1−h), p 2 (r 0 ) and p 2 (r) in the expression for G; each of these terms leads to a branch point, respectively at ω M , at k 0 corresponding to r + c (k 0 ) = r 0 , and at k r corresponding to r + c (k r ) = r.
Moreover, G possesses a pole at k 0 , which is exactly the non-modal pole referred to above, although there are no poles of G at ω 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 1 − h, and the pole and branch point at k 0 are not present when r 0 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 where R ± is the residue given in (3.11) evaluated using G ± , R + 0 (k 0 ) is the residue of the non-modal pole k 0 evaluated as if approached from above the branch cut, derived in appendix C.2 and given in (C 19) as the steepest descent integrals are defined as (3.15) and the jumps across each of the Steepest Descent branch cuts are calculated in appendix B to be 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.

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) which are amenable to asymptotic analysis in the limit x → ∞, 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 − iξ) ∼ Bξ ν + O(ξ ν+1 ) to leading order for small ξ with ν > −1, then Watson's Lemma implies that, for large x, where Γ is the Gamma function, and in particular, Γ (ν + 1) = ν! for integer ν. For each of the I q (x) integrals, this can then be interpreted as an algebraically decaying wave of phase velocity ω kq .
In order to find the decay rates of the steepest descent contours we are required to understand the behaviour of ψ 1 (r, k q − iξ) and ψ 2 (r, k q − iξ) for small ξ at r ∈ {1 − h, r, r 0 , 1}. Details of these can be found in appendix C. The result, given in equations (C 13) and (C 14), is that, for k = ω M − iξ, as ξ → 0 with ξ > 0, we find that and algebraically decaying like x − 5 2 when the source is within the region of uniform flow, and 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 ∆ G 0 as k → k 0 is derived in appendix C.2 as Finally, considering ∆ G r as k → k r , it is found in appendix C.3 that (3.20) By Watson's Lemma, this results in a wave convected with the flow speed U (r) and decaying algebraically like x −4 .
It may be noted that the decay rates for I 0 and I r are the same as predicted for a linear boundary layer flow profile by Brambley et al. (2012a). We now proceed to compare these results with previous literature. 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.

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)  The I 0 integral gives a wave with phase velocity equal to that of the mean flow at the location of the point mass source, U (r 0 ), provided the point mass source is in a region of sheared flow, r 0 > 1 − h. It can be observed in table 1 that agreement is seen in all three works for r 0 > 1−h. While Swinbanks did not consider the other cases in detail, this work finds further agreement for the I r contribution with Brambley et al.. In both the linear and quadratic shear flow cases, when the source is located within the region of sheared flow, the I 0 contribution is the slowest decaying term. When the source is located within the uniform flow region, the I ω M contribution is the slowest decaying term, although this is matched by I r contribution for linear shear. It should noted, however, that when r 0 > 1 − h we have in addition the contribution of the non-modal k 0 pole, which does not decay.
The difference in the behaviour of the I ω 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) In the case of linear shear flow, the U term is zero for k = ω M and the resulting expression is O(1) as k → ω M . In the case of a quadratic shear boundary layer, U is non zero and dominates A as k → ω M , providing a factor of (k− ω M ) − 1 2 . The remainder of the differences between the decay rates is explained, for I ω M , by the fact that (r in the quadratic shear case, where as for linear shear (r c − (1 − h)) ∼ (k − ω M ). For the I 0 and I r contributions, where we do not have r + c → 1 − h, and all the other terms are equivalent between the linear and quadratic cases, therefore giving the same eventual asymptotics scalings, although the pre-factors may vary significantly. Further details are given in appendix D.

Numerical Results
In this section, the above analysis is illustrated with some numerical examples. The Frobenius series solutions p 1 and 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 numerically expensive 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 p 1 and 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 |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 rotoralone tone at takeoff, while parameter set C might represent the same type of mode during the landing approach.
In section 4.1, we will explore the locations of the modal poles in the complex k-plane. This section will particularly focus on the k ± modal poles discussed in 3.3, including tracking these modal poles as they move behind the critical layer branch cut for certain parameters (by taking advantage of the ability of the Frobenius series solutions to analytically continue behind the critical layer branch cut). In section 4.2 we compare the various contributions from the critical layer branch cut described in section 3.3.3, including their large x behaviour, and show that these agree with the predicted large x behaviour from section 3.4. In section 4.3 the full solution in terms of x and r is plotted, and these results are compared to the linear boundary layer case. Finally, in section 4.4, we investigate how the results vary as we vary parameters, including the frequency ω, the boundary layer thickness h, the wall impedance Z, and the steady mean flow velocity M .

Pole Locations.
The locations of the poles in the complex k plane for parameter sets A1and A2are plotted in figure 5. In addition to the usual acoustic modes (denoted as * in figure 5), one k + and one k − pole is found for each parameter set. For parameter set A1, both the k + and k − poles are behind the critical layer branch cut, and so would not be found using conventional numerical methods, although the k + pole does still contribute to the total pressure field through the critical layer branch cut contribution, as described in section 3.3.3 above. In contrast, for parameter set A2, the k + pole is not behind the branch cut and takes the form of a standard modal pole, in this case a hydrodynamic instability surface wave. The stability of the modal poles is verified from the movement  Figure 5. Location of the poles in the complex k plane for parameter sets A1 (top) and A2 (bottom). Left: For real ω: acoustic modes with Re(k) < ω M ( * ); k + poles (+); k − poles (×); the critical layer branch cut (-); and branch points ω M and k0 for r0 = 1 − 9h 10 (•). Right: Trajectories of poles for −50 < Im(ω) < 0. Poles coloured red (left) and solid lines (right) denote poles contributing to the modal sum. Poles coloured blue (left) and dashed lines (right) denote poles hidden behind the branch cut (which varies with Im(ω)) and do not contribute.
of the poles in the k plane as Im(ω) is decreased from zero, following the Briggs-Bers Criterion (shown in the right-hand plots in figure 5); note that the critical layer branch cut also moves as a function of Im(ω). Of particular note is that the k + pole for parameter set A1 emerges from behind the critical layer branch cut as Im(ω) is reduced from zero, becoming a standard modal pole provided Im(ω) is taken sufficiently negative.
As discussed in section 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 in located above the branch cut and does not contribute towards the Fourier inversion. In the event that this k − pole were located below the branch cut, its contribution would almost exactly cancel the critical layer branch cut contribution, again seen by deforming onto the steepest 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 section 4.3 below.
Also plotted in figure 5 is the critical layer branch cut for k > ω M , together with the non-modal k 0 pole, which is only present for a point mass source within the boundary layer, r 0 > 1 − h. The effects of these non-modal contributions are illustrated in the next section. 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 on 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).

Branch cut contributions
We can further look 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 7(left)). Note that there is no non-modal k 0 pole when r 0 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 8 compares the numerically-computed steepest descent integrals with their predicted far-field rates of decay given in section 3.4, and a good agreement is seen in all cases.

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 (left) with the full solution including the critical layer and any unstable k + pole (right). When the k + pole is a convective instability, it clearly dominates the solution sufficiently far downstream, as it grows exponentially in x. In these near-field 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.
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.  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 r0 = 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. 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 only found to be unstable if the boundary layer is sufficiently thin. Even when both flow profiles give rise to a convective instability, we can see in figure 11 that the growth rate of the instability can be significantly different.
The change in nature of the k + pole in the quadratic case from convective instability to stable is clearly of significant importance. We therefore finally consider the variation of the solution as various parameters of interest are varied next. Fourier inversion, which also includes the k + pole. The parameter sets used from top to bottom are A1, B and C, with r0 = 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).

Variation of results with changing parameters
The variation of the acoustic modal sum is relatively well understood, so in this section we concentrate on the variation of the k + and k − modal poles as various parameters are varied. This includes whether Im(k + ) > 0, corresponding to a convective instability, or Im(k + ) < 0, corresponding to a stable modal pole hidden behind the branch cut that none-the-less contributes to the modal sum through the branch cut contribution. We also consider whether Im(k − ) > 0, meaning the pole is not included, or whether 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 Im(Z) being unstable and nearly hardwalled values of |Im(Z)| → ∞ being stable, as is seen from the k + poles movement in figure 12(b). The variation of stability as the frequency is varied remains unclear, although it appears likely from figure 12(c) that, for certain parameters, there would be a finite range of frequencies for which the k + pole would be unstable, while for frequencies either higher or lower than this range the k + pole would be stable. Note also from figure 12 that, in all cases, the k − pole is located above the branch cut and so does not contribute either to the modal sum or the critical layer branch cut.

Conclusions
In this work we have considered a cylindrical duct containing a parallel mean flow that is uniform everywhere except within a boundary layer of thickness h near the wall. Within the boundary layer, which need not be thin, the flow has a quadratic profile and satisfies the non-slip boundary condition at the duct wall, whilst maintaining a C 1 continuous flow. For such a flow profile, irrespective of the the thickness of the boundary layer, a solution to the Pridmore-Brown equation has been constructed making use of two Frobenius series expansions, valid for any wave number k. This enables the evaluation of the Greens function of the Pridmore-Brown equation, 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 Greens function solutions, which is equivalent to the solution for a point mass source in the linearised Euler equations. The Greens function is in some sense the general solution, as the solution subject to any forcing can be written as an integral over the Greens function, suitably weighted. Because of this, any behaviour the equations are capable of must necessarily be demonstrated in the Greens function solution, and so once the behaviour of the Greens function is understood, the equations cannot hold any further surprises. This is particularly important in this case, considering that the Greens function solution has been shown to include nonmodal contributions, such as the critical layer branch cut, which 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, where as in general the critical layer is due to a non-zero second derivative of the sheared flow profile. This also allows comparison to previous works, such as that of Swinbanks (1975) and (Félix & Pagneux 2007). Further, as the quadratic flow profile has a continuous first derivative, we have also been able to investigate the specific case of a point mass source at the boundary between uniform and sheared flow, r 0 = 1 − h, and we find this case retains the same behaviour as a point mass source within the region of uniform flow. In contrast, for the linear flow profile, r 0 → 1 − h is a singular limit. We therefore believe the results of the quadratic boundary layer flow profile to be in some ways typical of solutions for a general boundary layer profile.
The final solution for the Green's function for a point mass source was found to consist of a number of contributions. This solution is dominated, both upstream and, in the near field, downstream too, by the sum of modal poles. The modal poles, including acoustic and surface modes, are well known, and are typically used in mode-matching numerical methods. One complication found here to the surface modes is that a particular surface mode, here labelled k + , is found to sometimes disappear behind the critical layer branch cut (or, in other words, into the continuous spectrum). The contribution of this mode is not lost, however, and is in effect added to the critical layer contribution. In general, the modal poles present difficulty only in establishing which poles contribute upstream (x < 0) or downstream (x > 0) of the source. This can be established through application of the Briggs-Bers criterion, as summarised in section 3.3.
The effect of the critical layer, the focus of this work, always contributes downstream of the source, and is the dominant contribution to the far-field pressure downstream of the point mass source. This contribution, which results from integrating the Fourier inversion contour around the critical layer branch cut, may be viewed in three parts. The first contribution is from the k 0 non-modal pole, only present when r 0 > 1 − h, which does not decay with distance from the point source and is therefore dominant in the 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 steepestdescent non-oscillatory integrals I ω M , I r and I 0 , which are the results of accounting for the branch points coming from the critical points of the Pridmore-Brown equation. These contributions have a phase speed equal to, respectively, the uniform flow speed 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 ω.
For the linear flow profile boundary layer (Brambley et al. 2012a), a k − pole was found below and behind the critical layer branch cut that contributed to the critical layer. For the quadratic flow profile boundary layer here, this k − pole is always found to be above the critical layer branch cut, and so never contributes. We believe that this k − pole was an artifact of the unphysical linear boundary layer profile, although we have no direct way of demonstrating this. Incidentally, for the linear flow profile boundary layer, Brambley et al. (2012a) argued that the k + pole could never be behind the critical layer branch cut, as this would cause an unphysical discontinuity in the final solution; in fact, it is found here that when the k + pole is behind the branch cut, the unphysical discontinuity in the k + pole contribution is exactly cancelled by the I r steepest descent contour contribution, resulting in a continuous solution.
The various decay rates of the components of the critical layer have previously been predicted by Swinbanks (1975) and Brambley et al. (2012a); and a summary can be found in table 1. Swinbanks (1975) only considered the contribution from waves with phase velocity U (r 0 ), which are only present for a point mass source within the region of sheared flow, r 0 > 1−h. Swinbanks (1975) predicted these to behave as a constant amplitude plus a decay as O(x −3 ) in the far field. Brambley et al. (2012a) found the same result, despite Swinbanks considering a two dimensional flow in a rectilinear duct with an arbitrary flow profile and Brambley et al. considering only a constant-then-linear flow profile in a three dimensional cylindrical duct; in particular Swinbanks emphasised the importance of the non-zero second derivative of the mean flow profile, which is identically zero for a constant-then-linear flow profile. As a result, it would not have been unsurprising if these results were different. Here, the same result is again found, with the constant amplitude coming from the k 0 non-modal pole and the algebraic decay coming from the I 0 integral. This shows that this agreement is not by coincidence. For the critical layer contribution that propagates with phase velocity U (r) when r is within the boundary layer, we also find here the same result given by Brambley et al. (2012a) of an O(x −4 ) far-field decay.
The critical layer also contributes a term with phase velocity equal to the uniform flow velocity M , which is always present, and which dominates the other critical layer contributions in the far field whenever the point mass source is in the uniform flow region, r 0 < 1 − h. The amplitude of this term can decay at two different rates depending on the location of the source. When the source is within the uniform flow a decay rate of O(x − 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 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 accurately predicted using a surface wave dispersion relation (e.g. Brambley 2013), especially when the k + pole is located behind the critical layer branch cut; our experience in this work has been that it cannot, at least with the simplified surface wave dispersion relations that assume a thin boundary layer with a linear flow profile, although more complicated surface wave dispersion relations may prove more accurate. Another possibility for further investigation is to consider parameters on the stability boundary when the k + pole is neutrally stable. In this case, the k + pole is exactly located on the critical layer branch cut, and there would also exist a value of r 0 for which the non-modal k 0 pole and the k + pole overlap; this case has been explicitly excluded here. While this may seem a rather contrived case, a distributed sound source would correspond to an integral of source strength over all values of r 0 , and so k 0 and k + coinciding could be expected to occur for any parameters leading to exact neutral stability. One could also extend this problem to a non-constant mean density and sound speed making use of equation (2.5). For a given mean density profile ρ 0 (r) and sound speed profile c 0 (r), one could still construct a solution to the resulting Pridmore-Brown equation, taking careful notice of the potentially complex roots of c 0 (r). It would be possible to construct a solution using the Frobenius series solutions still so long as these are not double roots or of higher order, and ρ 0 (r) ρ(r) has at most regular singularities in the complex r plane. Except in the case where these singularities occur at r = r + c then the critical layer branch cut will still occur in identical form to that seen in our work, although the resulting scaling in the various limits seen above may vary. When retrieving these, 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).
Supplementary material. The Matlab source code used to produce the data plotted here is available at https://doi.org/10.1017/jfm.2022.753.
The coefficients a n and b n are then given by setting the remaining terms of the Laurent expansion of (A 4) to zero, resulting in the recurrence relation a n = 1 n(n + 3) k 2 a n−2 − k 2 M 2 h 4 a n−6 + 4Qa n−5 + 4Q 2 a n−4 − n−1 j=0 (−1) j (n + 2 − j)a n−1−j where we take a n = b n = 0 for n < 0. Requiring a 0 = b 0 = 1, we then find that and that b 3 is arbitrary, so we choose b 3 = 0. However, for the recurrence relation involving b 3 on the left hand side to hold, we also require that A is chosen to be (A 8) Here, the notation p c1 and p c2 denotes that these are two linearly independent solutions for p about the critical point r + c .
The Frobenius series solutions (A 6) are limited by a radius of convergence, in that the series converge if |r − r + c | < R for some radius of convergence R. Following from Fuchs Theorem (Teschl 2012, Theorem 4.5) this R is the distance between r + c and the next nearest singularity of the Pridmore-Brown equation, which is either at r = 0 or at r = r − c , and hence R = min |1 − h + Q|, 2|Q| .
(A 11b) with α n = β n = 0 for n < 0. These solutions are labelled p 11 and p 12 to indicate they are two linearly independent solutions to p expanded about the point r = 1. We now construct solutions to the homogeneous Pridmore-Brown equation p 1 (r) and p 2 (r) such that they are valid across the whole of [1 − h, 1]. We set p 1 (r) = p c1 (r) |r − r + c | < R A 1 p 11 (r) + B 1 p 12 (r) otherwise (A 12a) p 2 (r) = p c2 (r) |r − r + c | < R A 2 p 11 (r) + B 2 p 12 (r) otherwise (A 12b) First of all, note that these expansions are sufficient for a uniformly-valid expansion, as sketched in figure 14. Note also from figure 14 that the regions of convergence of the p c solutions and the p 1 solutions always overlap (except when k = ω M , which we exclude here). For any realr > Re(r + c ) contained within both regions of convergence, we may find the coefficients A 1 , A 2 , B 1 and B 2 are found by requiring continuity and continuous derivatives at r =r: A 1 A 2 B 1 B 2 = p 11 (r) p 12 (r) p 11 (r) p 12 (r) −1 p c1 (r) p c2 (r) p c1 (r) p c2 (r) (A 13) These coefficients A 1 , B 1 , A 2 and B 2 are independent of the specific choice ofr, and the resulting solutions p 1 and p 2 have not only C 1 continuity but C ∞ , since both are solutions to the Pridmore-Brown equation. In effect, p 1 analytically continues p c1 beyond its radius of convergence, and similarly p 2 analytically continues p c2 .
As described in (2.17), there is a jump in p c2 across the critical layer branch cut due to the log term in (A 6b). If the radius of convergence R is sufficiently large that r = 1 is within the radius of convergence, then no matching coefficients are needed, and this jump in p c2 obviously carries through to p 2 . In the other case, that R is sufficiently small that matching is needed, it follows thatr < 1. In this case, there is no jump in the matching coefficients A 1 , A 2 , B 1 and B 2 asr > Re(r + c ), and hence ∆ p 2 (r) = −2πiA p 1 (r)H(r + c − r).

(A 14)
This is analogous to the jump in p 2 given in (2.17), and shows that the jump in p c2 carries through the analytic continuation, as might have been expected a priori.
Whenř < 1 − h <r then the formula for ∆ G depends on whether ω M < k < k > or k > k > . In this case, we set ∆ G = ∆ G ω M for ω M k < k > , and ∆ G = ∆ G ω M + ∆ G > for k > k > , so that ∆ G > is the correction required for k > k > . In effect, ∆ G has two branch points, one at ω M and one at k > in this case, and by making this definition we may write By considering (B 1) in this case, we find that Finally, when we have 1−h <ř we must consider three cases: ω M < k < k < , k < < k < k > and k > < k. Similarly to the previous case, we consider ∆ G ω M = ∆ G for ω M < k < k < , and take ∆ G < and ∆ G > to be correction terms as k crosses k < and k > respectively. This leads to Appendix C. Asymptotic behaviours of G and ∆ G q In order to find the residue contribution of the non-modal pole at k = k 0 , and the decay rates of the steepest descent contours given by integrating ∆ G q along k = k q − iξ, we are required to understand the behaviour of p 1 and p 2 at r = 1 − h, r, r 0 , and 1 as k → ω M , k r , k 0 , where p 1 and p 2 are given in appendix A.3. Considering the evaluations at r, r 0 > 1 − h and 1, it can be noted that we must examine the cases that p 1 and p 2 are described as p c1 and p c2 respectively (as described in If instead r 0 > 1 − h, we find that Finally, setting k = ω M − iξ, so that Q = (1 − i)h M ξ/2ω + O ξ 3/2 (recalling that Re(Q) 0), we find that j(r 0 ) may be written to leading order as h 5 M 5/2 ω 3/2 1 − U (r 0 )/M (r 0 − 1 + h) 4 ξ 5/2 r 0 > 1 − h.
(C 14) C.2. Asymptotic behaviour as k → k 0 We now consider k → k 0 with r 0 > 1 − h. We have that Hence, in this limit, p 1 (r 0 ) and p 2 (r 0 ) may always be evaluated in terms of p c1 and p c2 , as we are always eventually within their radius of convergence. Hence, in this limit, For r = r 0 , the Bessel function, Hankel functions, and p 1 and p 2 all behave as O(1) quantities when evaluated at 1−h, r, and 1, resulting in O(1) behaviour for C 1 , D 1 , C 2 , D 2 , C 2 and D 2 . It can be shown that A = O(1) and that C.2.1. Behaviour of G as k → k 0 and the residue of the non-modal k 0 pole Substituting all the above into (3.3) (as k → k 0 from above) gives G(k) = −2M k 2 0 (ω − M k 0 ) 3πir 0 h 2 ω(k − k 0 ) 1 C + 1 D 2 − C 2 D 1 D 2 ψ 1 (r) r < r 0 D 1 ψ 2 (r) r > r 0 + O(1), (C 18) confirming a pole at k = k 0 that gives a residue contribution once integrated around of D 2 ψ 1 (r) r < r 0 D 1 ψ 2 (r) r > r 0 . (C 19) C.2.2. Behaviour of ∆ G 0 as k → k 0 Moreover, we may substitute all the above into ∆ G 0 from equations (3.16b) and (3.16c) to find the leading order contribution to ∆ G 0 as k → k 0 . First of all, we find the exact expression for ∆ G 0 to be (considering only r 0 > 1 − h, as otherwise ∆ G 0 ≡ 0) A p 1 (r 0 ) C − 1 D 2 − C 2 D 1 + 2iπAD 1 D 2 × D 2 ψ − 1 (r) r 0 > r D 1 ψ − 2 (r) r 0 < r.
(C 20) Using asymptotics above, to leading order we find that (C 21)