## Abstract

The downstream propagation of high-frequency acoustic waves from a point source in a subsonic jet obeying Lilley's equation is well known to be organized around the so-called ‘cone of silence’, a fold catastrophe across which the amplitude may be modelled uniformly using Airy functions. Here we show that acoustic waves not only unexpectedly propagate upstream, but also are organized at constant distance from the point source around a cusp catastrophe with amplitude modelled locally by the Pearcey function. Furthermore, the cone of silence is revealed to be a cross-section of a swallowtail catastrophe. One consequence of these discoveries is that the peak acoustic field upstream is not only structurally stable but also at a similar level to the known downstream field. The fine structure of the upstream cusp is blurred out by distributions of symmetric acoustic sources, but peak upstream acoustic beaming persists when asymmetries are introduced, from either arrays of discrete point sources or perturbed continuum ring source distributions. These results may pose interesting questions for future novel jet-aircraft engine designs where asymmetric source distributions arise.

## 1. Introduction

Over the past 60 years, much research into the acoustic pressure field arising from jet flow has been based on the acoustic analogy originally developed by Lighthill [1,2] and later extended by Lilley [3]. The Lilley analogy models the acoustic field due to fine-scale turbulence in a subsonic jet by propagating the acoustic field of equivalent multipole-type convected compact sources through a mean parallel shear flow [4]. This allows a close approximation of the acoustic field from slowly spreading subsonic jets, facilitating the reduction of the governing propagation equations to a single third-order linear wave equation—Lilley's equation.

In the subsonic regime, the linearity of the Lilley equation may be exploited by constructing a Green function which describes the flow–acoustic interaction effects from a simple monopole source. The total sound field at a particular frequency is then determined by convolving the Green function with a specified source distribution [4].

Even in the flow geometries where the Green function is separable, its numerical determination is non-trivial. However, where the wavelength of propagation from a compact acoustic source is much smaller than the scales of the mean flow (as is the case within a few diameters of a round jet nozzle where turbulence is at its peak), a high-frequency approach may facilitate the calculation. Although the approximate crossover to this regime will differ depending on the flow, such solutions have been shown for subsonic single-stream parallel jets to agree with numerical computation down to Strouhal numbers of 1 by Tester & Morfey [5] and more recently as low as 1/2 by Wundrow & Khavaran [6].

In the case of high-frequency propagation, there are two popular asymptotic approaches. First, when the flow is radially symmetric an approximation of the Green function by separation of variables leads to a series of truncated Fourier modes each containing a leading-order Wentzel–Kramers–Brillouin (WKB) expansion in the large wavenumber [4,6,7]. A first term approximation of this sum forms the well-known MGB method due to Mani *et al.* [8]. However, abandoning non-zeroth mode terms oversimplifies the acoustics by losing information regarding the directional behaviour of jet noise. In this paper, we show that retaining these terms is vital to an accurate calculation of the directivity of the Green function three-dimensional free space.

Alternatively, the Green function may be described using ray theory or geometrical acoustics [9,10]. Recent developments have shown that one can now seek full-scale ray solutions to the Green function field in a generic flow [11,12], including the important effects of complex rays [13]. The trade-off between a computationally efficient symmetric WKB expression and a full-asymmetric solution is not automatically required [11,12]. In those papers, we provide a practical numerical algorithm to calculate acoustic fields in the presence of generic subsonic flows (symmetric or asymmetric), using both real and complex rays. We exploit that machinery here to give a glimpse of the rich underpinning mathematical structures of the acoustic field from a point source.

Much of the previous work for subsonic jets using the linearized Lilley equation considers the Green function in the far-field downstream of a point source. The downstream region of the jet contains the acoustic ‘cone of silence’ [4,14], the shadow zone of the high-frequency Green function. The field within the cone of silence is generated from complex rays only. It is delineated from the rest of the acoustic field by a region of peak amplitude.

Near to the loci of the peak acoustic field, the ray approximations generate a singular, divergent amplitude. This is due to the tangency of pairs of rays emitted from the point source [15] causing the vanishing of a Jacobian [11,12] in the denominator of the local approximation. The effective boundary of the cone of silence is thus a fold caustic of the rays, being the simplest member of the catastrophe hierarchy of bifurcation sets ([15], [16], ch. 36, [17,18]). It is well known that the field across this fold catastrophe may be uniformly approximated using Airy functions [6,15].

In this paper, we initially consider the far-field Green function of a point source in a typical subsonic jet flow in a high-frequency (high-Strouhal number) limit using ray approximations in all three-dimensional directions. Importantly, we include both upstream and downstream analysis. We restrict our analysis to subsonic flows governed by ideal and perfect gas laws and do not consider instability or shock-associated noise that requires a partial or full nonlinear treatment.

We find that the downstream fold catastrophe is complemented by a region of peak noise extending forward from an off-axis point source in the flow, in an upstream direction. The intersection of the associated caustic of rays with the surface of a far-field sphere centred on the point source is a cusp. The cusp is formed by three interacting real rays and is the next most complicated caustic in the hierarchy of catastrophe theory. We show that a uniform approximation of the Green function generates an upstream field amplitude comparable in magnitude to that of the downstream fold catastrophe. For a point source, the cusp is decorated by rapid oscillations of the Green function caused by interferences of the three real rays, uniformly approximated locally by Pearcey functions [16].

While upstream beaming of supersonic screech and broadband noise is well established [19,20], this involves fundamentally different mechanisms from those discussed in this paper. For subsonic flows emanating from round jets with symmetric source distributions, upstream beaming is non-existent [21] and the downstream sound pressure levels dominate angular spectra. However, if an asymmetry is present in the source distribution then upstream beaming may appear. Some indications of the upstream beaming we observe here may be found in the theoretical and experimental results of Powles & Tester [22], wherein the shielding effects of a subsonic jet were studied by artificially generating a monopole noise source external to the flow. However, the upstream angular range they cover and the experimental sampling frequency within that range are limited. Hence they do not cover the entire domain where the full complexity of the upstream structures that we have calculated here are revealed.

It should be emphasized that we study a fully three-dimensional system in the absence of boundaries. The present structure contrasts with previous appearances of cusps in aeroacoustic problems; for example, where caustics are focused downstream in subsonic flows and formed by rays reflecting off two-dimensional boundaries [23]; where cusps are embedded in hyperbolic umbilic catastrophes inside ducts caused by edge diffraction [24], or in supersonic flows by the focusing of two-dimensional supersonic shock waves [25]; and to model downstream shock cells in jets [26].

Although the downstream acoustic field has been well studied [6,14,27], our analysis illustrates additional novel structures within that regime. The fold catastrophe itself is actually embedded within a swallowtail catastrophe, the appearance of which is controlled by the radial source distance from the jet centreline. This is the third most complicated structure in the catastrophe hierarchy [16]. The local form is of further mathematical interest as it arises from the interaction of four rays, but unusually two of which are complex.

As members of the catastrophe theory hierarchy these caustics are structurally stable, i.e. resistant to perturbations. They can, therefore, be assumed to be present for a generic range of parameters. We show elsewhere that the cusp persists for a range of subsonic Mach numbers [12].

However, this is not the complete story, as any realistic acoustic field modelled by the linear Lilley analogy may be approximated as the sum of a distribution of such point sources [4]. For asymmetric discrete and continuous distributions we show that, although the fine structure oscillations in the Green function are washed out, the catastrophe structures generating local peaks of field amplitude still visibly persist in the field. Such washing is also likely to occur in other realistic scenarios.

The outline of the paper is as follows. In §2, we revisit the high-frequency modal solution of Wundrow & Khavaran [6], the standard numerical testbed for point source off-axis jet noise. Extending the modal solution to account for more complicated turning-point structures, we plot a global solution (§3) to reveal the catastrophe structure using a prototypical axisymmetric jet profile. In §4, we begin our analysis of the cusp catastrophe in the upstream direction. Ray trajectories are visualized and classified according to the Poisson summation method of Wundrow & Khavaran [6]. Analysis then moves on to the downstream region in §5, whereby a local diffraction integral is derived to describe the cone of silence boundary, resolving ray conservation issues found in [6,11]. The final sections of this paper consider the impact of both upstream and downstream singularities, starting with a direct magnitude and energy calculation in §6, and culminating with a source distribution study based on correlation integrals in §7.

## 2. Green function solution

In this section, we present the modal Green function solution to Lilley's equation, based on an extension to Wundrow & Kharavan [6]. Although this solution relies on separable flow it nevertheless acts as a benchmark for our more general ray analysis.

The parallel shear flow configuration specifies an infinitely long cylindrical jet according to the flow profiles
*M*(*r*) and *a*(*r*) are the acoustic Mach number and sound speed ratio, respectively. All profiles are functions of the radial coordinate *x*,*y*,*z*) of figure 1, so that *x*-axis with **i** the unit vector in that direction. Variables (2.1) take on ambient values

As in [6] for computational convenience we consider two Green function solutions of Lilley's equation given by
*ω* is defined as
*x*,*r*,*φ*), and *x*_{R}=(*x*,*r*,Δ*φ*) and point source *x*_{s}=(*x*_{s},*r*_{s},*φ*_{s}) is shown in figure 1,
*φ*_{s}=0 for single-source computations without loss of generality.

The solution proceeds by determining *G*_{ω} first, then obtaining the ultimate solution *x* and periodic nature of the field in Δ*φ* makes a combined Fourier transform series in axial wavenumber *k*_{1} and azimuthal index *n* possible
*G*_{n} satisfy the reduced equation
*Φ*=(1+*κM*(*r*))/*a*(*r*), *a*_{s}=*a*(*r*_{s}), *κ*=*k*_{1}/*k*_{0} and *v*_{j} (*j*=1,2) is a solution to the homogeneous equation
*r*_{>} branch of (2.8) subject to the conditions
*k*_{1} of (2.6) in the far-field regime, *k*_{1}=*k*_{1★}=*k*_{0}*κ*_{★} [29]. As shown in [6], this Fourier inversion leads to the solution of

The computation of *v*_{1} and *v*_{2}. Numerical routes are possible; however, WKB solutions offer both a closed-form estimation and our first glimpses at the inner workings of the point source field.

### (a) WKB solutions

Solutions of (2.9) can be formed using the WKB method [7,30], neglecting *v*_{1,2} as linear combinations
*r*_{δk} is a reference turning point, one of a number *T* real turning points, which satisfy *Q*_{n}(*r*_{δk})=0, ordered according to *r*_{δT}>*r*_{δT−1}>⋯>*r*_{δ1}≥0. These turning points partition the positive real line into subdomains. The solution in each subdomain is glued to the next using the connection constants *k*) denotes the solution to the left of the *r*_{δk}th turning point). Turning points are caustics of the modes, and their existence and position should not be confused with the ray caustics we analyse in later sections. From (2.10), we see that turning points vary according to the flow (as *q* depends on *Φ* the function that characterizes the flow), modal index *n*, wavenumber and receiver position *θ* (via the stationary point *Φ*, the techniques which we use to deal with the turning points will be the same.

For most *θ*∈[0,*π*), *n* in commonly used smooth subsonic flow profiles *a*(*r*) and *M*(*r*), there exists only one turning point (*T*=1) on the positive real line [8]. For example, in figure 2 we show the root loci for the isothermal flow used in this paper: *M*(*r*)=*M*_{J} sech^{2}(2*r*), *M*_{J}=0.9, *a*(*r*)=1. By drawing a line of *θ*=const. up to a particular modal curve, there is typically only one intersection, i.e. one turning point. Consequently, the most useful form is the simple turning-point problem, and for a large range of angles (*θ*<8*π*/9 (160°), *θ*≠*θ*_{S} for *n*=0, for the flow used here) one can use the results of [6], where the solutions of *v*_{1},*v*_{2} (*m*_{rδ}=1 (see (A 8) and appendix A). This determines the Wronksian

The use of Ai on its own in (2.16) implies that it is ‘transitionally’ uniform at *r*_{δ1}, i.e. it does not lead to a divergence in the modal field, but is the first term of the uniform expansion theory applicable to ordinary differential equations derived in appendix A. Numerically, we find these are appropriate for *modal* calculations even in the event of higher order turning points.

In the event of higher-order and multiple turning points, we must amend (2.16) only, as *v*_{1} is never non-uniform at *r*_{δ1}. These complicated turning-point geometries are not accounted for in the literature as they are typically found in upstream regions of the jet hitherto neglected. The exception to this is the higher-order turning point in the *n*=0 mode for downstream *θ* discussed below. This occurs in a variety of subsonic jet flows [4–6] and is key to understanding a new phenomenon presented in §5.

First, we deal with higher-order turning points that occur in (2.16). For the flow regime used in this paper, figure 2 shows that there are two significant angular regions, *θ*=*θ*_{S} for mode *n*=0, and *θ*≥8*π*/9 for modes 1≤*n*≤14, where the turning points become second order. For the angle *θ*_{S} this occurs due to a coalescence between a real positive turning point and its image on the negative real axis, and those for *θ*>8*π*/9 occur due to a complex conjugate pair *r*_{δ2,3} coalescing on the real line.

The second-order turning-point structure at *θ*_{S} for *n*=0 is easy to elicit from the local form of *q*^{2}(*r*,*θ*) about the turning point at *r*=0, noting that *q*^{2} is even in *r* for this flow regime and truncating at second order gives
*θ*_{S} yields *T*_{0}(*θ*_{S})=0, so that *θ*_{S} as the ‘Schwarzschild line’. The second-order turning point exists in (2.20) when *n*=0.

In both regions *θ*=*θ*_{S} and *θ*≥8*π*/9, we can use the local uniform form (A 4) from appendix A with second-order condition *m*_{rδ}=2 to determine the transitional solution. However, for the modal problem, as opposed to ray analysis in §5, we use the alternative parabolic cylinder function representation of the relevant Bessel function ([16], appendix A and ch. 10.16)

For large angles *θ*>8*π*/9, multiple turning points exist as a result of the complex conjugate pair of turning points, *r*_{δ2,3}, scattering along the real axis post coalescence and introducing a new region of decay (unphysical growing solutions are discarded) between *r*_{δ1}<*r*′<*r*_{δ2}. This is shown in figure 2 as *θ*=const. lines intersect the modal curves three times. In this case rather than derive a new solution we prefer to modify the simple turning-point solution using the *θ*=*π* as the infinite number of turning points—corresponding to trapped rays [12]—lies beyond our modification.

## 3. Combined catastrophes

In this section, we display the magnitude of the modal solution, *φ*∈[0,2*π*) and *θ*∈[0,*π*) for the first time. This reveals key information concerning the structures about which the whole field is organized. Figure 3 shows the modal solution for an isothermal jet, *a*(*r*)=1, using flow profile *M*(*r*)=*M*_{J} sech^{2}(2*r*), *M*_{J}=0.9, source point *x*_{s}=(0,0.75,0), and far-field radius *R*=200. The field is expressed in terms of both Strouhal number *St*=*k*_{0}*r*_{J}/*πa*(0)*M*_{J} with *r*_{J}=1/2, *a*(0)=1, and wavenumber *k*_{0}. As in [6], a modal criterion of *n*_{max}=*O*(*k*_{0}*r*_{s}) is used to truncate (2.13), replacing limits *n*=±*n*_{max}.

The modal field (2.12) shown in figure 3 for this particular subsonic flow is typical of off-axis source positions, *r*_{s}>0. Figure 3 also exhibits the benefits of a two-dimensional plot with extended plot range as opposed to constant Δ*φ* cross-sections shown in [6]. It is now clear from the diffraction patterns alone that

We can locate the caustic structures that organize these diffraction patterns by using the ray tracing apparatus of [11,12]. The caustics, where the ray amplitude is singular, can be located by calculating the loci of points where the ray Jacobians vanish. The local catastrophe forms can be deduced by counting the number of singular rays involved at the caustic.

As figure 3 shows there are two caustics inhabiting two distinct regions of the field, a fact reinforced by the spherical polar plots of figure 4. The first in the downstream region, *θ*<*π*/2, organizes the diffraction pattern delimiting the cone of silence—the region of exponential attenuation existing for low polar angles (see also [11]). To date this has been classified as a fold caustic, controlled and decorated locally by the Airy function. This last fact is easy to see when we correct for ray non-uniformity using the local form
*A*_{CE} and *c*_{0,0} are analogous to the ray phase and amplitude, and critically *ξ*_{1} controls the essentially one-dimensional oscillations. The local form (3.1) is analogous to, but distinct from, (2.16) used to provide uniformity in the modes.

Under the ray approximation the scenario is simple: two real/complex rays coalesce together (or are tangential in physical space) violently at the caustic giving birth to a complex/real ray pair in a manner equivalent to the bifurcations of two polynomial roots. However, as we shall show in §5 of this paper this caustic is more than just a fold. There are two peculiarities suggesting otherwise: first, the existence of a caustic evaporations point where the caustic appears to disappear from the field; second, apparent cusp-like kinks just visible in figure 3 above *θ*=*θ*_{S} near to Δ*φ*=*π*/2,3*π*/2. As we show in §5, these kinks are dependent on source parameters.

The second caustic of *θ*>*π*/2, and sees us reap the rewards of retaining non-zeroth modal terms in (2.13). The cusp, a new discovery of *P* [[16], eqn (36.2.14)]. We can show this as before by evaluating the uniform expansion for the case of three rays
*c*_{0,0} and *A*_{CE} are as before—though take on different numerical values—and the presence of two control variables *ξ*_{1,2} implies the oscillations vary along two principal axes. As is well known, the cusp unfolds into two fold caustics, mimicking the roots of the polynomial *t*^{3}+*ξ*_{2}*t*+*ξ*_{1} where three rays/roots coalesce on *ξ*_{1}=*ξ*_{2}=0, and two coalesce on the bifurcation set *θ*=0 by mapping the exact solution found in ([[16]], ch. 36.5).

In both cases, the caustics shown in figure 3 are calculated at fixed radius, *R*. If we allow *R* to vary it follows that the fold caustic forms a smooth sheet (modulo the cusp kinks) and the cusp forms a cusp rib, both propagating into the near field.

While both caustics form the key structural components of the field, there are two additional structures that are important. The first is the Schwarzschild line *θ*_{S}, characterized by a double turning point in the radial equation which is independent of Δ*φ* (see (2.21) and figure 2). Not only do fold oscillations wash out upon crossing *θ*_{S} (figure 3), but for the flow studied here the caustic evaporation marks the tangential intersection of the fold with *θ*_{S}.

The second structure, the ray delimiter shown by the white line in figure 3, is not really a structural entity *per se*, but a useful tool for our local analysis in §5. This structure delimits the region where the continuation of the ray pair that coalesce at the fold is ‘indirect’. By this we mean that both ray trajectories initially propagate away from the receiver before being turned by the flow [11,32]. The key to this region, used in §5, is that it contains the caustic evaporation points and, as we discover, additional cusps.

In the following sections, we shall analyse each of the structures independently before assessing their influence upon the field in both isolated point source cases and as contributors to a source distribution.

## 4. The cusp and its far-field genesis

We start by analysing the structure of the cusp. Here our point of view concerns the genesis of the cusp structure formed by those rays that propagate out to the far field as opposed to the analysis of intricate near-field structures that arise from a trapped ray field [12]. Analysis is aided by use of the Poisson sum formulation of [6] shown in appendix B, recasting the modal summation (2.13) as *ν*_{★}, equivalent to rays [13], so we can write ray contributions as the sum
*ν*_{★}, and ray trajectories [9] to label the rays using the Poisson index, *m*_{★}. This provides closure to the large *θ* ray contributions missing from [6], and describes the rays’ behaviour from an analytical perspective.

The ray trajectories are numerically integrated solutions of the initial value problem in Cartesian coordinates (*x*,*y*,*z*) as shown in [11,12]
** x** and

**are the ray position and slowness vector (or wavefront normal), respectively, both functions of travel time**

*p**τ*with initial conditions

**(0)=**

*x*

*x*_{s}and

**(0)=**

*p*

*p*_{s}. Additionally,

**is the unit matrix, ⊗ denotes an outer product and ∇**

*I*_{x}is the gradient operator in (

*x*,

*y*,

*z*).

As in [11,12], (4.2), (4.3) are recast into a two-point boundary-value problem allowing us to cut the singularity structure of the cusp and record the relevant trajectories. As we know from the unfolding of the cusp, we need only cut it in two places to observe two and three ray caustics, the latter of which will confirm we really are dealing with a cusp. Cutting the cusp is achieved by parametrizing the receiver *x*_{R} by *θ*, with Δ*φ*=const. generating cross-sectional slices of

The results of our two-point boundary calculation are shown in figures 5 and 6 with accompanying *ν*_{★} bifurcation diagrams for off-symmetry line and symmetry line cases, respectively. In the first of these figures, we make several observations. The first is that the coalescing ray couple share the same orientation of path about the jet axis, while the third ray takes a path encircling the axis in the opposite sense. The choice of orientations is reflected in the Poisson index, *m*_{★}, with coalescing rays sharing the same value.

The second observation concerns the qualitative behaviour of the real ray triplet appearing at angles *θ*>*θ*_{c}, where *θ*_{c} is the cusp intersection angle. The majority of these rays, including those that are singular on the cusp, initially propagate in the direction of the flow (figure 5) before creeping back against the flow. This counterintuitive behaviour is one of the reasons the cusp has not been undiscovered until now. Analytically this upstream propagation can be seen as the result of a battle between flow and slowness vector ** p** in (4.2). Initially, the flow term

**dominates before conceding to the term**

*M***⋅**

*T***as the ray propagates. This means that**

*p***, which maintains a similar value to**

*p*

*p*_{s}the initial orientation of the ray trajectory, ultimately prevails.

The final observation relates to the orientation of ray trajectories about the jet axis when the receiver is placed on the symmetry line, as shown in figure 6. As opposed to the off-axis case there is no preferred route for the coalescing rays: two rays propagate with opposite orientations about the axis; the third propagates through the axis. This discord between real ray trajectories is reflected in the Poisson index for *θ*>*θ*_{c} (figure 6*c*).

## 5. Ray vanishing and apparent fold evaporation

### (a) A ray vanishes

Having previously considered the cusp we now look at the fold caustic near *θ*_{S}. As figure 3 shows, the fold appears to evaporate near to *θ*_{S}, where its oscillations also fade away. The field for *θ*<*θ*_{S} is supported by two real rays that coalesce at the fold, where for *θ*>*θ*_{S} one of these real rays has either vanished or become complex. A local analysis confirms the latter.

Before we go into detail, we can summarize what is known about the Schwarzschild line. First, as it is a radial equation property, it is independent of Δ*φ*. Second, on this line it is known that an indirect saddle, that is, a saddle *ν*_{★} of the integral (setting *m*=0 in (B 2))
*θ*_{S} represent a locus of second-order turning points as *ν*=*ν*_{★}→0 in (2.20) (since *ν*=*n*/*k*_{0}), but ray evaluations of (5.1) erroneously vanish as *θ*→*θ*_{S} according to

It can be seen by the ray delimiters in figure 3 that the fold caustic rays are both saddles of (5.1) in the vicinity of the evaporation point, meaning that we only need to analyse the integral *r*_{δ}=*r*_{δ1} from now on.

### (b) Local analysis

In this section, we show that the ray solutions about the Schwarzschild line and caustic evaporation point are controlled by the sum of a quadratic polynomial and logarithmic term. This form can be seen as the unfolding of the local form valid on the Schwarzschild line as the receiver is perturbed away from it.

We start by reconsidering the radial equation on the Schwarzschild line noting that its locus is characterized by the saddle *ν*_{★}=0 and a second-order turning point at *r*=0 (figure 2). The local behaviour of *q*^{2}(*r*,*θ*) is given by (2.18), (2.19), so we can write the local radial equation as
*ν*=*n*/*k*_{0} of appendix B. Now under the ray approximation (4.1) we know that the radial solutions must be evaluated at the saddle point, i.e. *ν*=*ν*_{★}=0 in (5.3), due to the recasting of (2.13) using the Poisson formula. As before, evaluation on the Schwarzschild line gives *Q*^{2}(*r*|*θ*_{S})∼*T*_{2}(*θ*)*r*^{2}. In this case, (5.3) is soluble in terms of Bessel functions of order ±1/4, (see (A 4)). For example, for *v*_{2}(*r*_{s}) we have
*v*_{1}(*r*) as *q*^{2} tends to a constant in the far field and the above approximation would violate this. The solution is to change the phase reference (i.e. the domains of integration) of *v*_{1,2} by examining their joint phase contributions to the integrand of *ζ*(*r*_{s} | *ν*)+*ζ*(*r* | *ν*). We partition these contributions by rewriting the limits of integration as *r*_{δ}<*r*′<*r*_{s}, and take the out-of-flow approximation, *M*(*r*)→0, for integration over

We reflect this change of integration by recasting our radial Green function components as

Crucially, the function *r*_{s} and *r*_{δ} interactions only, whereas *r*_{δ}. Under the same Taylor approximations as before we can solve *ν*=*ν*_{★}=0 in terms of Bessel functions. However, rather than rewriting (5.4) again, we use the special case of the Pearcey function on its symmetry line ([16], eqn (36.2.19)) to give the local radial structure on *θ*_{S} as

The local form (5.7) suggests the strategy we now use to understand the underlying structure of the wavefield. We further our local analysis by including non-zero *ν* terms in an expansion that unfolds the Schwarzschild line behaviour. Again we use the Taylor expansion (2.18) in *M*→0 in *ν*^{2}/*r*^{2}. The phase integrals are now
*θ*_{S}, where the turning point is *r*_{δ}=0, the turning point of (5.9) satisfies the closed-form expression
*Q*(*r*_{δ} | *ν*)=0 under the truncated Taylor regime.

With both phases (5.8), (5.9) integrable, we expand their solutions sequentially using the small parameters *T*_{0},*ν*→0, as *θ*→*θ*_{S}. We perform the expansion in *T*_{0} to first order and then expand in *ν* up to powers of *ν*^{4}, in what is a post-paraxial approximation [33] and one motivated by the Pearcey function in the radial solution (5.7).

Once expanded the phase of *P*_{h} of (5.2) can be written
*f*(*r*_{s},*θ*,*M*_{J}) collecting terms independent of *ν* and

Retaining only those terms that affect the dynamics of the local form (i.e. only functions of *ν*) leads us to discard *f*(*r*_{s},*θ*,*M*_{J}) from (5.11). We can now define a new codimension-3 integral form based on the expansion (5.11), which takes on a similar guise to those of catastrophe theory once rescaled to leave *ν*^{4} independent of receiver variables. Modulo the constants due to the Jacobian and logarithm of this rescaling, and the slow function (*i*/2*πr*_{s}*Q*(*r*_{s}|*ν*))^{1/2} of (5.1), we have
*φ*,*r*_{s},*θ* for a constant Mach number, *M*_{J}=const.

The new integral form allows us to make a straightforward numerical comparison straight away. Initially, we choose *r*_{s}=0.45 to ensure that the underlying approximation (5.3) makes a decent estimate of the phases *ζ*, where the structure of *r*_{s}=0.75, pertinent to figure 3, is considered below. Figure 7 shows both modelled (5.13) and modal (2.13) fields for *r*_{s}=0.45 with caustic and Schwarzschild structures. In the case of the modelled field, we note that the caustics are centred about Δ*φ*=*π*/2, the zero of *I*.

The first comparison we make is upon the qualitative features of the field. As expected from local expansion about *ν*=0, the best match is found in the vicinity of the evaporation point where both fold saddles tend to 0. The agreement is notably broken by the increasing misalignment of the model's oscillations with increasing lateral distance from Δ*φ*=*π*/2. There are two reasons for this: first, the local mapping becomes increasingly poor as the saddles move away from *ν*=0; second, the function *I* is only the transitional solution about the evaporation point and a full mapping would include derivatives of *I* with respect to

The second comparison we make is of the structural behaviour. We note that the model predicts a continuous branch of caustics, caustics 1 and 2 in figure 7, whereas only caustic 1 can be determined in the modal solution using the methods of [11,12]. As we show below, caustic 2 is resolved by examining the Riemann sheet structure of *I*. The curvature of caustic 1 and that of the modal field match near to the evaporation points, while the Schwarzschild lines match exactly. Note that for this source radius neither field contains the cusp kinks of figure 3. We now go on to explain when and why they occur.

The appearance of caustics exhibited by figure 7 can be understood by examining the underlying singularity structure of *I*. This is given in *P*_{h}/∂*ν*=∂^{2}*P*_{h}/∂*ν*^{2}=0. These conditions are equivalent to two fourth-order polynomials, and, as we show in Stone *et al.* [35], yield a singularity structure that is the swallowtail catastrophe ([16], ch. 36.4) albeit with complex coefficients. Here the caustic is given implicitly by
*canonical* and *structurally stable*. This may be confirmed by a formal coordinate transformation [34].

The singularity set (5.14) is shown in figure 8*a* with equivalent cross-section for *r*_{s}=0.45. The swallowtail structure shows how the fold caustics, 1 and 2, of the model appear and how they merge with the Schwarzschild plane *θ*_{S}. An important feature of figure 8*a* is the existence of cusp ribs for *b* for swallowtail cross-sections at *r*_{s}=0.75,0.9, in the (Δ*φ*,*θ*) plane. This leads us to conclude that the cusp-like kinks near to the evaporation points of figure 3 really are cusps. A notable prediction that is true even though *r*_{s}≥0.75 is perhaps outside the range of validity of the original Taylor expansion (2.18).

Even though the model field has a symmetric caustic structure about the evaporation point, this is clearly not reflected in the field amplitude, as it would be by the swallowtail catastrophe [[16], ch. 36.3]. The intricacies of the model field and its implications for *I* as *π*/2,*θ*_{S}). We do this for *r*_{s}=0.45, though the ray mechanisms behind the Schwarzschild line and evaporation point are the same for all source radii considered in this paper. We examine further the case *r*_{s}=0.75 in [35].

Using the control points ** a**={

*a*

_{1},…,

*a*

_{7}} of figure 7, we determine which of the saddles contribute to the field using the method of steepest descents (SD) [29]. For convenience we unfold the logarithmic singularities of (5.13) by transforming

*ν*=

*e*

^{w}, where the saddles

*w*

_{★}satisfy the saddle condition

This transformation returns the primary saddles *w*_{★} and an infinite number of images, i.e. *w*_{★,k}=*w*_{★}+2*πik*, *k*=0,±1,±2,…. Note that *k*=0 returns the primary saddle. The SD contour then ascends from the valley at complex infinity with orientation chosen in order to ensure convergence. We label the saddles (ℓ,*k*) where ℓ labels one of the four primary saddles, with *k* as before. The contributing saddles are then mapped back to *ν*_{★} and *I* is evaluated using the ray sum
*ν*_{★}→0 as *θ*→*θ*_{S}, in agreement with the order of *ν*_{★} found in the vanishing ray amplitude of

Figure 9 shows the SD paths for control points {*a*_{1},…,*a*_{6}}. Starting with *a*_{1} in figure 9*a*, we note that the receiver, located in the cone of silence, picks up a complex ray (2,0) and an infinite number of exponentially small images of itself and of its conjugate (4,0). The conjugate pair then coalesce at *a*_{2}, coincident with caustic 1 (the fold, *θ*_{f}; see figure 9*b*), as do all their image pairs. As the receiver moves from *a*_{2} to *a*_{3} the complex rays emerge post coalescence with (2,0),(4,0) a real ray pair. So far, so canonical.

Control point *a*_{4} (figure 9*d*) sees all saddles (4,*k*), ∀*k*, move off and cluster at infinity. This is of course equivalent to *ν*_{★}→0 at *θ*_{S}. Moving to *a*_{5} then sees the return of (4,*k*) to the finite plane, making exponentially small contributions in the presence of (2,0), which remains real. Note that going from *a*_{4} to *a*_{5} the ray expansion has undergone a Stokes phenomenon picking up exponentially small ray contribution (3,0).

Traversing around the evaporation point to *a*_{6} lying on caustic 2, the ray expansion again diverges. This time the caustic involves exponentially small rays (3,ℓ+1),(1,ℓ),ℓ≥0, which lie on the branch cut of *I* in the *ν*-plane (which we take from *I*. While the caustic structure is that of the swallowtail, half of this structure pertains to caustics of exponentially damped rays. This caustic would appear in

We can examine the evaporation point using *a*_{7} with SD paths shown in figure 10. Whereas ordinarily crossing *θ*_{S} one ray and its images make excursions to infinity, the evaporation point sees two rays, (1,0) and (4,0), coalesce while making the same excursions. This results in an otherwise canonical caustic being damped at infinity in the presence of real ray (2,0).

We can sum ray contributions in *ν*_{★} using figure 11*a* for Δ*φ*=1.3, extending the bifurcation diagram of [6]. Recording the new complex branches, e.g. (4,−1) the dominant complex ray for *θ*>*θ*_{S}, we plot the ray field of *I* in figure 11. This figure shows how the rays erroneously predict a kink at *θ*_{S} due to vanishing amplitude, where the uniform correction is *I* itself.

While *I* models the modal field locally about the evaporation point and confirms a swallowtail structure, the number of additional rays present in the upstream region means that

## 6. Cuspoid magnitude and energy analysis

Now we have determined fully the controlling structures of the Green function we are in a position to examine the relative importance of each. The first comparison we make is a direct amplitude and energy calculation. In both cases, we consider the point source solution.

The aim of this comparison is to identify the largest contributions of the field as a function of Δ*φ*. As is well known, the uniform asymptotic field about a canonical caustic is greatest not on the caustic but in the region just behind it. In light of this, figure 12 shows the modal field, *φ*=*π*, we traverse the Schwarzschild line.

It can be seen that, while the swallowtail cross-sections are larger for the most part in Δ*φ*, the maximum values associated with each caustic are commensurate, and though the cusp drops to half the value of the fold away from the maximum, it has a larger distribution in Δ*φ*. This suggests that we should take into account the distribution of the caustics and their attendant diffraction wavefield, not just their maximum values. A suitable way of doing this is using an energy analysis: integrating the energy flow at constant radius *R* over a representative region containing the caustic maximum values.

The energy of the Green function is proportional to the amplitude squared, so we will consider normalized (by *R*^{2}) surface integrals of the form
*D* is the region of integration in (Δ*φ*,*θ*), and gives the total energy, *E*_{T}, when *D*=[0,2*π*)×[0,*π*). In the following, we use for *D*: disc domains *B*_{i}={** p**∈(Δ

*φ*,

*θ*):|

**−**

*p***|≤1} about the evaporation point (**

*q**i*=1,

**=(**

*q**π*/2,

*θ*

_{S})), and cusp point (

*i*=2,

**=(**

*q**π*,

*θ*

_{c})); and rectangular domains

*R*

_{1}=[0,2

*π*)×[0,

*π*/2),

*R*

_{2}=[0,2

*π*)×[

*π*/2,

*π*), for downstream and upstream fields, respectively.

The results of this study are shown in figure 13, where we have included energy due to *B**=*B*_{1}∪*B*_{2}, which is equal to the energy sum of discs about both evaporation points (counting *B*_{1} twice) and the cusp point. These results show that the energies contained in equivalent area discs about the cusp and the evaporation point are of similar magnitude, and that the combination of these energies using *B** accounts for over half the field energy. This is more than that of both the upstream and downstream regions individually.

The relevance of energy contained in *B*_{i} is due to the beaming of the caustics towards observers in aircraft noise evaluations. For instance, when the source is placed at *φ*_{s}=*π*/2, *B*_{2} is a useful measure of the cusp's impact at flyover certification points [36], despite the cusp rib in three dimensions providing a more complicated acoustic signature at a microphone than analysis at *R*=const. As we show in §7, distributions of sources lead to distributed caustics, meaning they are relevant at all microphone positions, particularly lateral and flyover positions where the jet engine is at, or is near to, its maximum power setting.

## 7. Aggregate effect of singularities

To test the persistence of the caustic structures, in this section we consider the effect of multiple off-axis sources upon the catastrophe fields analysed above (using the same flow profile). We first consider the combination of discrete source distributions as one might expect in engineering implementations of acoustic analogies to non-separable flows [12]. In such circumstances the convolution of the Green function with a source distribution is replaced by a series of weighted point source solutions of the fields similar to that shown in figure 3. We then move on to variations of continuum ring source models popular in aeroacoustics, as used in [5,6].

### (a) Multiple point sources

The effect of multiple point sources on the field can be examined easily if they lie on the same source radius due to the separable nature of the Green function (2.6). A sum of *N*_{s} Green functions with equal source weighting is given by
*φ*=*φ*−*φ*_{s1}, Δ*φ*_{1m}=*φ*_{s1}−*φ*_{sm} and *φ*_{sm} is the *m*th of *N*_{s} source azimuths as illustrated in figure 14*a*(i).

The aggregate fields for *N*_{s}=2,3 with radius *r*_{s}=0.75 are shown in figure 15. Unsurprisingly, in the cases where these sources are symmetrically placed, i.e. figure 15*a*,*c*, the fields are periodic in Δ*φ*. Conversely, asymmetric source positioning in figure 15*b*,*d* leads to an asymmetric field which is somewhat disordered. In each case, the field is complicated by the overlapping diffraction patterns, effectively leading to a completely decorated field. One notable consequence of this distribution is the apparent reduction in cusp intensity. We might then ask if continuing the process of adding point sources in the manner of (7.1) *ad infinitum* (

### (b) Ring sources

Ring sources are frequently used in aeroacoustics as a model simplification (e.g. [5,37]). They reflect the turbulence physics by concentrating the source in the shear layer of the jet. We use the correlation integral setting of [12] to arrive at our ring source directivity function *V* ′,*V* ′′ are integrations over real space, ** η** are separation coordinates, that is, the vector displacement between two source points,

**=**

*η*

*x*_{s}′′−

*x*_{s}′, say.

The ring source model proceeds from (7.2) by taking the source correlation as
*M*_{ij} are source position terms so that *x*_{⊥s}′=(*y*_{s}′,*z*_{s}′), and *δ*(*x*_{s}′−*x*_{s}) reflects the far-field and *x*-independent result of the source distribution. The correlation effects are reflected in *δ*(** η**), which exhibits zero correlation between any two source points that are not coincident. This can be seen as the limit of Gaussian source models used in jet noise [38,39]. We make one further simplification, and that is we limit ourselves to just one component of the inner product of the integrand in (7.2), i.e. we set

*M*

_{ij}=

*M*(

*x*_{⊥s}′), ∀

*i*,

*j*. This makes the source independent of frequency. We claim that the conclusions we draw about the impact of singularities upon

#### (i) Symmetric ring source

We start by examining a perfectly symmetric ring source. Taking *M*(*x*_{⊥s}′)=*δ*(*r*_{s}′−*r*_{s})/*r*_{s}′ we have
*δ*_{nm} is the Kronecker delta). This result equals that of [6] and is shown in figure 16*a*. Not only is it invariant in Δ*φ*, but it also exhibits the effects of repeatedly overlapping the cusp and swallowtail caustics giving a constant boundary to the cone of silence and completely washing out the cusp. The latter observation answers the question posed in §7a, as the cusp effects are eliminated in the

#### (ii) Weighted ring source

The washing out of the cusp shows that perfect symmetry of the ring source eliminates any asymmetry present in the Green function. However, it remains to be seen how the cusp impacts the field if there are small perturbations to this symmetry. In this section, we consider perturbations to the ring source weight.

The weighted source uses source correlation of the form *M*(*x*_{⊥s}′)=*h*(*φ*_{s}′;*φ*_{P},*α*_{T})*δ*(*r*_{s}′−*r*_{s})/*r*_{s}′ allowing for the presence of azimuthal variability and its tuning through parameters *φ*_{P} and *α*_{T}. The former centres the asymmetry, while the latter concentrates it. This results in an integral of the form

Choosing a Gaussian form for the weight, i.e. *a*(ii) for an example of weighting). We then find the azimuthal component is given by

An example of *b* for various *α*_{T}. The effect of increasing asymmetry, i.e. increasing *α*_{T}, affects *n*=0. Conversely, *α*_{T}→0 allows

Equation (7.7) is computed for *φ*_{P}=*π*/2 and *α*_{T}=5 as shown in figure 16*b*. Representing quite a large deviation from the symmetric ring source field, figure 16*b* draws a stark comparison with the symmetric field of figure 16*a*. Not only do we see the re-emergence of the cusp upstream, but the total field resembles a Green function with smoothed diffraction patterns.

#### (iii) Perturbed ring source

In the case of an asymmetric source position, the argument of the delta function can be adapted. The resulting integrals, however, are more cumbersome than before. Some progress can be made if a perturbed ring source is considered. This permits the construction of a perturbation expansion in a small parameter, *ϵ*≪1, say, as shown in figure 14*a*(iii). As a function of *ϵ* the delta function then expands as

The directivity factor *O*(*ϵ*), which can be written as

The effects of the perturbations on the ring source are shown in figure 17. This figure shows that even for small changes in the asymmetry of the source there are significant changes in the qualitative behaviour of the field. The upstream beaming is immediately recognizable as soon as the field is perturbed, with the field smoothed as in the case of the weighted ring sources. However, unlike the weighted source cases, the largest cusp-related magnitudes are of the same order as the swallowtail-related diffraction patterns delimiting the cone of silence. The overall impact of these overlapping swallowtail singularities is reduced with increasing *α*_{T}, with the concentration of the diffraction field again resembling a point source field in its organization.

Given the studies undertaken in this section it is not difficult to surmise that coupling both weighted and perturbation cases together would probably see the cusp as non-negligible.

## 8. Conclusion

The discovery of the upstream acoustic cusp in a parallel subsonic shear flow with asymmetric source distributions provides a new perspective on how high-frequency acoustic waves propagate in jets. This fully three-dimensional effect leads to a distinctive diffraction pattern (in the absence of any scattering from structures) with significant peak amplitude in the upstream direction. For a point source of acoustic radiation, the cusp behaves canonically, being decorated locally by the Pearcey function, and is formed by an interacting ray triplet creeping back against the flow.

Spurred on by this discovery we revisited the downstream field organized about a fold caustic and Schwarzschild line, the latter characterized by a radial potential of the form used in black hole theory. We resolved the interaction of these structures by performing a local analysis that introduced a new controlling form, not simply Airy as with folds, but a complex Pearcey-like function. Although appearing non-canonical, the overall controlling structure is that of a swallowtail catastrophe and is thus locally structurally stable. Not only does this new form predict cusps downstream arising from swallowtail cross-sections, explaining the kinks in the fold at large source radii, but we are also able to predict new complex ray shedding as one traverses the Schwarzschild line, with the effect of smoothing the field.

With the cusp and swallowtail structures resolved, it is a simple question of which forms dominate when the acoustic field is the product of multiple sources. While equally weighted point sources suggest the total elimination of the cusp, perturbations to a ring source, representative of the ‘imperfect’ round jet, introduce asymmetries restoring the upstream beaming effect.

Currently, there is a lack of experimental data comprehensively covering the upstream angular regions for asymmetric source distributions in subsonic flows that form the subject of this study. As trends in aircraft design suggest closer coupling between the engine and the aircraft structures, leading to an effective asymmetry in effective acoustic source distributions, a systematic experimental study is highly desirable. If the presence of upstream caustics were to be observed, this might have practical implementations for novel jet–wing or jet–nozzle acoustic interactions in future aircraft design.

## Data accessibility

Computational documents (CDFs) of Green function and ray trajectories are available within the electronic supplementary material. Additional data supporting this paper are openly available from the University of Southampton repository at http://dx.doi.org/10.5258/SOTON/403450 and http://dx.doi.org/10.5258/SOTON/403414.

## Competing interests

We have no competing interests.

## Authors' contributions

All authors participated in the formulation of the mathematical model and wrote the paper. J.T.S. did the numerical simulations. All authors gave final approval for publication.

## Funding

This work was supported by EPSRC (EP/M508147/1), Rolls-Royce and Mathematical Sciences (University of Southampton).

## Acknowledgements

The authors thank the reviewers for their valuable comments and suggestions. J.T.S. thanks Rolls-Royce for its support.

## Appendix A. Extensions to simple turning-point problem

**(a) Transitional and uniform forms**

In this section, we outline the theory used to derive the ‘transitional’ uniform solutions valid in the vicinity of an isolated turning point *r*_{δ} (where *r*_{δ} is any *r*_{δk}; see §2) . By ‘transitional’ we mean the leading-order term with error of magnitude

Starting with the radial equation
*r*_{δ}, of order *m*_{rδ} if *v*, but is less recondite. Making the transformations *r*=*r*(*η*) and *d*/*dη*), giving
*V* satisfies
*V* has two linearly independent (transitional) solutions *V* _{±} in terms of Bessel functions *J*_{±1/(2+mrδ)}
*A*_{0}=const. determined from (2.11). We then have *V* given by (A 5).

The Bessel functions (A 4) can be expressed in terms of Airy functions when *m*_{rδ}=1 and parabolic cylinder functions when *m*_{rδ}=2 ([16], ch. 10.16). In both cases, it is expedient to absorb *k*_{0} and (−1)^{mrδ} into *η* of (A 3) to give

The uniform forms derived here are equally applicable to radial equations evaluated at *s*=*r*_{s} when the WKB solutions (2.14) are singular, i.e. when *r*_{s}=*r*_{δk} of §2.

**(b) Three turning points**

In the event of three turning point, *T*=3, we modify the expansion of (2.16), *T*=1 problem. Note that, in all cases, we can solve for *v*_{1}(*r*) by making the replacement *r*_{δ1}→*r*_{δ3} in (2.15).

The first two subcases are *r*_{δ3}<*r*_{s} and *r*_{δ2}<*r*_{s}<*r*_{δ3}; we simply re-use (2.16) with *r*_{δ1}→*r*_{δ3}, ignoring *r*_{δ1,2}. However, in the subcases *r*_{s}<*r*_{δ2} we must take care to integrate the phases over the regions of decay 0<*r*′<*r*_{δ1} and *r*_{δ2}<*r*′<*r*_{δ3}. To demonstrate we use the most common subcase that occurs in figure 2 for the sources used in this paper, i.e. *r*_{s}≥0.45, in which case *r*_{δ1}<*r*_{s}<*r*_{δ2}. Each phase must acquire a decaying factor from integration over *r*_{δ2}<*r*′<*r*_{δ3} and this is easily accounted for in phase *p*.1 of (A 9) by replacing *r*_{δ1}→*r*_{δ2}, and in phase *p*.2 by integrating from *r*_{s} to *r*_{δ1} and then from *r*_{δ1} to *r*_{δ2}. This leads to
*r*_{s}<*r*_{δ1}.

## Appendix B. Poisson summation and saddle point asymptotics

The Poisson summation introduced in [6] recasts (2.12) by decomposing (2.13) according to
*r*_{s}>*r*_{δ}=*r*_{δ1}, for the radial solutions. This gives
*ν*=*n*/*k*_{0} and function *α*_{m}=Δ*φ*+2*πm* of the Poisson index *m* have been introduced, along with notational changes *ζ*_{n}(*s*)→*ζ*(*s*|*ν*), *Q*_{n}(*s*)→*Q*(*s*|*ν*), where *s* is defined in (2.14). We use

The saddle point expansion of (B 2) is given as
*ν*_{★} are saddle points satisfying the saddle point condition

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3770105.

- Received December 1, 2016.
- Accepted April 18, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.