1. Introduction
In many applications such as rocket engines, imperfectly expanded supersonic jets are often accompanied by a powerful emission of tonal sounds. Such tones are commonly referred to as jet screech. In addition to the intensified noise emission, it may also lead to disastrous structural damage because of sonic fatigue. It is, therefore, of practical importance to understand the mechanism of supersonic jet screech and to devise effective ways to suppress these screech tones.
Screech occurs in imperfectly expanded supersonic jets, and such jets are often characterized by quasi-periodic shock cells and are complex in nature. Jet screech was first discovered in an experiment conducted by Powell in the 1950s. In his pioneering work, Powell (Reference Powell1953a) proposed the well-established feedback loop which consists of four stages, i.e. the instability growth in jets, the interactions between shock and instability waves, the acoustic waves propagating upstream and the receptivity of the shear layer at the nozzle lip. Powell proposed the phase and gain conditions which must be satisfied to sustain the feedback loop. For the phase condition, the frequency of the fluctuation 
$f$ is supposed to close the feedback loop as
\begin{equation} \frac{N}{f}=\frac{d}{U_{c}}+\frac{d}{c^{*}}+\varPsi, \end{equation}
where 
$d$ denotes the distance between the nozzle lip and the sound source, 
$U_{c}$ is the average velocity of the instability waves travelling downstream, 
$c^{*}$ is the speed of sound propagating upstream, 
$\varPsi$ represents an additional phase delay and 
$N$ is an integer. For the gain condition, the gain from each of the four stages must satisfy
\begin{equation} Q\eta_{s}\eta_{u}\eta_{r}\geq1, \end{equation}
where 
$Q$ denotes the gain associated with the growth of the instability waves and 
$\eta _{s}$, 
$\eta _{u}$ and 
$\eta _{r}$ represent the efficiencies of energy transmission in the last three stages, respectively. The resonance conditions were then reconsidered by examining the energy exchange between the instability and acoustic waves at the sound source location and the nozzle exit (Landau & Lifshitz Reference Landau and Lifshitz1958). In recent work, these conditions (Landau & Lifshitz Reference Landau and Lifshitz1958) have been rewritten in terms of magnitude and phase conditions, details of which can be found in Jordan et al. (Reference Jordan, Jaunet, Towne, Cavalieri, Colonius, Schmidt and Agarwal2018) and Mancinelli et al. (Reference Mancinelli, Jaunet, Jordan and Towne2021).
Two important characteristics of jet screech have been widely studied over the past few decades. The first is the screech frequency. Considering the feedback enhancement during the loop, Powell proposed that the screech frequency 
$f$ can be calculated via
\begin{equation} f=\frac{U_{c}}{s(1+M_{c})}, \end{equation}
where 
$s$, 
$U_{c}$ and 
$M_{c}$ are the shock cell spacing, the convection velocity and the convective Mach number of the instability waves, respectively. Subsequently, Tam, Seiner & Yu (Reference Tam, Seiner and Yu1986) proposed the weakest link theory suggesting the screech as the limit of the broadband shock-associated noise when the observer angle approaches 
$180^\circ$. In 1999, Panda (Reference Panda1999) discussed the link between screech and hydrodynamic-acoustic standing waves, and a new formula was developed. In early measurements (Powell Reference Powell1953a), it was found that the screech frequency experienced abrupt changes as the inlet pressure increased. This frequency jumping phenomenon is commonly referred to as mode staging, and four different stages, i.e. stages A, B, C and D, were observed by Powell, among which stage A can be further divided into two stages named 
$\rm {A}_{1}$ and 
$\rm {A}_{2}$ (Merle Reference Merle1956). It was found that the azimuthal mode of both sound and instability waves changes as mode staging occurs, which shows a strong sensibility to the facility and initial conditions (Anufriev et al. Reference Anufriev, Komarov, Kuptsov, Mel'nikov and Sergienko1969; Gutmark, Schadow & Bicker Reference Gutmark, Schadow and Bicker1990; Panda, Raman & Zaman Reference Panda, Raman and Zaman1997), and the switch from one mode to another is nearly immediate (Nagel, Denham & Papathanasiou Reference Nagel, Denham and Papathanasiou1983). Despite different stages showing different characteristics, it is interesting to note that they can appear simultaneously in one jet flow (Raman Reference Raman1996). To interpret mode staging, Shen & Tam (Reference Shen and Tam2002) suggested that it was the neutral acoustic waves, rather than free acoustic waves, that complete the feedback loop in 
$\rm {A}_{2}$ and B modes. Recent works (Edgington-Mitchell et al. Reference Edgington-Mitchell, Jaunet, Jordan, Towne, Soria and Honnery2018; Gojon, Bogey & Mihaescu Reference Gojon, Bogey and Mihaescu2018; Li et al. Reference Li, Zhang, Hao and He2020; Mancinelli et al. Reference Mancinelli, Jaunet, Jordan and Towne2021) showed that both the 
$\rm {A}_{1}$ and 
$\rm {A}_{2}$ modes were closed by the neutral acoustic waves (or guided-jet modes). Gao & Li (Reference Gao and Li2010) used the original phase condition shown in (1.1) and inferred the value of 
$N$ in their numerical study. They showed that 
$N$ differed across various stages, and when this difference was considered the prediction was in very good agreement with experimental data. However, the mechanism behind this mode transition and its high sensitivity to the initial conditions are yet to be clarified. In addition, it was found that nonlinearity could arise during mode staging in circular jets (Mancinelli et al. Reference Mancinelli, Jaunet, Jordan and Towne2019) and screech tones were not independent but were instead nonlinearly phased locked to each other in rectangular jets (Walker & Thomas Reference Walker and Thomas1997), which increases difficulties for modelling.
In addition to the screech frequency and mode staging, noise directivity is the second characteristic that has been widely studied. It was found that acoustic radiation at the fundamental frequency appeared strongest in the upstream direction, whereas at the harmonic frequency there was a strong beaming to the side of the jet (Powell Reference Powell1953a). To explain this, Powell proposed the monopole array theory, which was generally in good agreement with experiment results. The directivity of the fundamental tone and its harmonics in supersonic round nozzles was measured by Norum (Reference Norum1983), the results of which were compared with the monopole array theory when nine monopoles with a parabolic intensity distribution were considered. One particularly interesting observation was that the strongest emission appeared somewhere near 
$150^\circ$, not 
$180^\circ$, to the downstream jet axis. A quick decay occurred when the observer angle approached 
$180^\circ$, which could not be predicted by Powell's model. As a matter of fact, if all the sound sources were monopoles and the frequency of the fundamental tone was obtained by (1.3), the acoustic radiation would be the strongest at 
$180^\circ$. Following Norum's idea, the directivity pattern of three equal-spaced monopoles of various intensities was studied by Kandula (Reference Kandula2008). It seemed that this variation did not affect the location of the main directivity lobe. The influence of temperature on the directivity pattern (Massey et al. Reference Massey, Ahuja, Jones III and Tam1994) was also investigated in round jets, but no significant difference was found from the unheated one. In the case of rectangular nozzles, the screech problem may become more complicated compared with axisymmetric round jets. However, this problem can be greatly simplified when the rectangular nozzle is of high aspect ratio, in which case it reduces to a two-dimensional problem. Numerical simulations (Berland, Bogey & Bailly Reference Berland, Bogey and Bailly2007) and experiments (Walker & Thomas Reference Walker and Thomas1997) were conducted to study the directivity patterns of the screech and its harmonics in rectangular jets of high aspect ratios and the results were similar to those of Powell (Reference Powell1953a). Recently, Tam & Parrish (Reference Tam and Parrish2014) considered another nonlinear interaction mechanism between shock, instability and acoustic waves. They then proposed a model to predict the lobe position in the directivity patterns. The result was in good agreement with experimental data (Norum Reference Norum1983) at harmonic frequencies, but appeared less so for the fundamental tone, in particular for the lobe position.
As argued by Powell, four stages were involved in the screech cycle. Among the four stages, it is believed that the interaction between the shock and instability waves plays a critical role in understanding the physics of screech. Not only because this interaction produces sound that is directly measurable, but also because it is the key to understanding the noise generation mechanism. Despite its importance, not many theoretical models have been proposed to predict the interaction. The reason is in part due to the complex flow nature present in the interaction (Manning & Lele Reference Manning and Lele2000), especially when the shock waves are intense in highly underexpanded and overexpanded jet flows. Harper-Bourne & Fisher (Reference Harper-Bourne and Fisher1973) used Powell's phased-array model to study the interaction between the disturbance in jet shear layers and shock cells. The frequency of the emitting sound was obtained. Subsequently, Tam (Reference Tam1987) developed a shock cell model composed of time-independent waveguide modes (Tam, Jackson & Seiner Reference Tam, Jackson and Seiner1985). The turbulence structures were modelled using a noise initiated at the shear layer at the nozzle lip. The weak interaction between these two components then gave rise to the sound field. Although the sound generated by shock–vortex interaction was analytically studied in the first part of the work, it was remarked by the author that (Tam Reference Tam1987) ‘enormous amounts of numerical computations are required’ for practical calculations. Thus, ‘a model source function’ was used instead in light of the extreme complexity of the practical evaluation of the formulation to calculate the directivity patterns of broadband shock-associated noise. Lele (Reference Lele2005) further developed Tam's theory. Based on the method proposed by Lighthill (Reference Lighthill1952), he used the wave packet model to describe the instability waves initialized by white noise, and a vortex-sheet model was utilized to calculate the shock cell structure. These two components were inserted into the Helmholz equation as the source term. The sound field was obtained by integration. This model was used in subsequent numerical simulations (Wong, Jordan & Honnery Reference Wong, Jordan and Honnery2019). However, the sound sources in the above-mentioned models were obtained by a simple combination of the shock and instability waves. A correct source term directly from governing equations would be more desirable. A somewhat different approach to model the sound generation is the so-called shock leakage mechanism. It was proposed by Manning & Lele (Reference Manning and Lele1998) and Manning & Lele (Reference Manning and Lele2000), and theoretically developed by Suzuki & Lele (Reference Suzuki and Lele2003) and Shariff & Manning (Reference Shariff and Manning2013). Recently, it was experimentally observed by Edgington-Mitchell et al. (Reference Edgington-Mitchell, Weightman, Lock, Kirby, Nair, Soria and Honnery2021b). In addition, a very recent work (Nogueira et al. Reference Nogueira, Jordan, Jaunet, Cavalieri, Towne and Edgington-Mitchell2022b) numerically studied the linear stability characteristics of shock-containing jet flows, where the shock was assumed to be of small amplitude and a sinusoidal form. It was found that the characteristics of the instability waves in shock-containing jets are different from those in shock-free supersonic jets, and a new interpretation of screech was proposed based on this observation.
In 1994, Kerschen & Cain (Reference Kerschen and Cain1995) developed an analytical shock–instability wave interaction model for two-dimensional planar vortex-sheet flows. The source term was obtained directly from the governing equations. One shock cell was considered to interact with the instability waves near the vortex sheet. However, it was found that the radiation field peaked at 
$48^\circ$ to the downstream jet axis, which contradicted experimental observations. Despite numerous attempts, an analytical and quantitative study of the interaction between shock and instability waves, which is capable of predicting not only screech frequencies, but also directivity patterns of screech tones, is yet to be seen. This paper aims to develop such a model to predict the sound arising from the interaction between shock and instability waves. The model follows the asymptotic expansion method proposed by Kerschen & Cain (Reference Kerschen and Cain1995), but a more realistic jet and shock cell structures are considered.
This paper is structured as follows. Section 2 presents a detailed analytical derivation of the model, while § 3 shows the prediction of the screech frequency and the directivity patterns of the fundamental tone and its harmonics. The near-field pressure and noise generation mechanism are subsequently discussed. Conclusions are presented in § 4.
2. Analytical formulation
2.1. The interaction model
To enable analytical progression, we start with a vortex-sheet model. As shown in figure 1, the coordinate axes 
$(x^\prime,y^\prime )$ are chosen to be parallel and perpendicular to the nozzle centreline, respectively. Here 
$D$ is the jet height (note that 
$D$ is generally not equal to the height of the nozzle, and we take the height of the fully expanded jet as the base flow height (Tam Reference Tam1972)). Velocity 
$U_{1}$ is the jet velocity at the nozzle exit plane while 
$U$ is the velocity of the fully expanded jet after exiting from the nozzle. The fully expanded base flow described in figure 1 takes the form
\begin{equation} \boldsymbol{u}_{0}= \begin{cases} 0, & |y^\prime|>D/2, \\ U \boldsymbol{e}_{x^\prime}, & |y^\prime|\leq D/2, \end{cases} \end{equation}
where 
$\boldsymbol {e}_{x^\prime }$ is the unit vector in the 
$x^\prime$ direction. We assume that the shock and instability waves are of small amplitudes and can be linearized around the base flow and described by linear theories. Of course, it is known that the interaction between the shock and instability waves primarily occurs several shock cells downstream from the nozzle exit (Suda, Manning & Kaji Reference Suda, Manning and Kaji1993; Kaji & Nishijima Reference Kaji and Nishijima1996; Malla & Gutmark Reference Malla and Gutmark2017). At these locations, instability waves are likely to grow to a significant amplitude where nonlinear effects become important and the instability waves may start to saturate and even decay. However, it is known that linear theories can predict the wavelength of these large coherent structures well beyond the linear stage (Crow & Champagne Reference Crow and Champagne1971; Jordan & Colonius Reference Jordan and Colonius2013; Edgington-Mitchell et al. Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021a). We may, therefore, use linear stability analysis to determine the wavelength (hence convection velocity) of instability waves, which are particularly important for the generation of screech tones. The linear growth rate calculated describes the early evolving of the instability wave and its role is discussed separately in subsequent modelling. The interactions between the shock and instability waves several shock cells downstream are likely to be nonlinear in strongly underexpanded jets. However, a linear interaction model may suffice to describe the interaction between weak shock and instability waves. Besides, similar to the successful prediction of the wavelength and convection velocities of large coherent structures by linear theories, a linear theory may still possess the many essential features of a nonlinear jet screech. We therefore start with a linear interaction between the shock and instability waves. With these assumptions we start to seek an analytical model describing noise generation due to the interaction between shock and instability waves.
Figure 1. Schematic of the vortex-sheet flow configuration and Cartesian coordinates. The origin is fixed at the centre of the nozzle while 
$x^\prime$ and 
$y^\prime$ represent the streamwise and cross-flow coordinates, respectively.
Given the fact that the Reynolds number is high and the instability waves are essentially inviscid, we start from the Euler equations shown as follows:
$$\begin{gather} \frac{\mathrm{D} \rho}{\mathrm{D} t^\prime } + \rho\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$
$$\begin{gather}\rho \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D} t^\prime}={-}\boldsymbol{\nabla} p, \end{gather}$$
$$\begin{gather}\frac{\mathrm{D} s}{\mathrm{D} t^\prime} = 0, \end{gather}$$
where 
$t^{\prime }$ denotes time, 
$\boldsymbol {u}=(u,v)$ the velocity, 
$p$ the pressure, 
$\rho$ the density, 
$s$ the entropy and 
$\mathrm {D}/\mathrm {D} t^\prime =\partial /\partial t^\prime + \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }$. Because the entropy increase across a weak shock is a high-order small term (Kerschen & Cain Reference Kerschen and Cain1995), the isentropic condition is used here.
To determine the solution, both kinematic and dynamic boundary conditions need to be satisfied across the vortex sheet. The dynamic boundary conditions read
\begin{equation} \left.p_{+}\right|_{y^\prime = h^{'}} = \left.p_{-}\right|_{y^\prime = h^\prime}, \end{equation}while the kinematic boundary condition requires
$$\begin{gather} \left. v_{+}\right|_{y^\prime = h^\prime} =\left(\frac{\partial h^\prime}{\partial t^\prime}+ \left.u_{+}\right|_{y^\prime = h^\prime}\frac{\partial h^\prime}{\partial x^\prime}\right), \end{gather}$$
$$\begin{gather}\left. v_{-}\right|_{y^\prime = h^\prime} =\left(\frac{\partial h^\prime}{\partial t^\prime}+ \left. u_{-}\right|_{y^\prime = h^\prime}\frac{\partial h^\prime}{\partial x^\prime}\right), \end{gather}$$
where 
$({\cdot })_+$ and 
$({\cdot })_-$ represent the quantities outside and within the jet flow, respectively, and 
$h^\prime$ denotes the displaced 
$y^\prime$ of the vortex sheet.
Following Kerschen & Cain (Reference Kerschen and Cain1995), we use 
$\delta$ and 
$\epsilon$ to denote the strength of the shock and instability waves, respectively. These two parameters are assumed to be of the small magnitude representing small perturbation compared with the mean jet flow. Note that although shock waves can be intense in strongly imperfectly expanded jets, they could also be infinitesimally weak when the wave angle approaches the Mach angle (Anderson Reference Anderson2017). In this model, we consider a slightly imperfectly expanded jet, in which case the shock-associated perturbations can be relatively weak, and thus a linear model may be developed (Tam & Tanna Reference Tam and Tanna1982; Tam et al. Reference Tam, Jackson and Seiner1985; Kerschen & Cain Reference Kerschen and Cain1995). The velocity field can be expanded using these two parameters as (Kerschen & Cain Reference Kerschen and Cain1995)
\begin{equation} \boldsymbol{u}=\boldsymbol{u}_{0}+\delta \boldsymbol{u}_{m} +\epsilon\boldsymbol{u}_{v}+\delta^{2}\boldsymbol{u}_{m2} +\epsilon^{2}\boldsymbol{u}_{v2}+\delta\epsilon\boldsymbol{u}_{i}+\cdots, \end{equation}
where 
$\boldsymbol {u}_{0}$ is the mean velocity, 
$\boldsymbol {u}_{m}$ represents the linear perturbation due to shock waves and 
$\boldsymbol {u}_{v}$ is the linear unsteady perturbation due to instability waves. For higher orders, the 
$\delta ^2$ term represents nonlinear steady modification of the shock waves and is independent of time. As we are interested in sound generation, this term can be neglected. The second-order term 
$\epsilon ^2$ represents the nonlinear correction to the linear stability waves, and we neglect it by only considering the leading-order contribution from the 
$\epsilon$ term. The omission of this term follows immediately if 
$1 \gg \delta \gg \epsilon$ is assumed, which implies 
$\epsilon ^2$ is the highest-order term of the expansion. The 
$\delta \epsilon$ term represents the interaction due to shock and stability waves, as a result of which sound is generated. Physically, this entails that the shock and instability waves that interact to produce sound in realistic flows may be approximated by the linear shock and instability solutions to the base flow. The pressure, density and vortex-sheet displacement have similar expansions, i.e.
$$\begin{gather} p=p_{0}+ \delta p_{m} +\epsilon p_{v}+\delta^{2}p_{m2}+\epsilon^{2}p_{v2}+\delta\epsilon p_{i}+\cdots, \end{gather}$$
$$\begin{gather}\rho=\rho_{0}+ \delta\rho_{m}+\epsilon\rho_{v}+\delta^{2}\rho_{m2}+ \epsilon^{2}\rho_{v2}+\delta\epsilon\rho_{i}+\cdots, \end{gather}$$
$$\begin{gather}h^\prime=h^\prime_{0}+ \delta h^\prime_{m} +\epsilon h^\prime_{v}+\delta^{2}h^\prime_{m2}+\epsilon^{2}h^\prime_{v2}+\delta\epsilon h^\prime_{i}+\cdots. \end{gather}$$
The mean temperature across the jet flow is different, while the mean pressure 
$p_{0}$ remains identical in both regions. If a perfect gas is assumed, then 
$c_0=\sqrt {\gamma p_0/\rho _0}$ can be used to calculate the mean speed of sound inside and outside the jet, and it is straightforward to show that 
$\rho _{0-}c_{0-}^{2}=\rho _{0+}c_{0+}^{2}$. We define 
$M_{-}=U/c_{0-}$ and 
$M_{+}=U/c_{0+}$ to denote the mean Mach numbers inside and outside the jet, respectively. Substituting all the expansions to the Euler equation and boundary conditions, and collecting the terms 
$O(\delta )$, 
$O(\epsilon )$ and 
$O(\delta \epsilon )$, we obtain the equations governing the shock and instability waves and their interaction, respectively.
2.2. The shock model
The 
$O(\delta )$ terms in the governing equations representing the linear perturbation induced by the shock wave satisfy
$$\begin{gather} \frac{\mathrm{D}_{0}p_{m}}{\mathrm{D} t^{'}}+\rho_{0}c_{0}^{2} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{m}=0, \end{gather}$$
$$\begin{gather}\rho_{0}\frac{\mathrm{D}_{0}\boldsymbol{u}_{m}}{\mathrm{D} t^{'}}={-}\boldsymbol{\nabla} p_{m}, \end{gather}$$
where 
$\mathrm {D_0}/\mathrm {D} t^\prime$ denotes 
$\partial /\partial t^\prime + \boldsymbol {u}_0\boldsymbol {\cdot }\boldsymbol {\nabla }$. These two equations can be combined to yield
\begin{equation} \boldsymbol{\nabla}^{2}p_{m}-M^{2}\frac{\partial^{2}p_{m}}{\partial x^{\prime2}}=0. \end{equation}
Equation (2.14) reduces to the Laplace equation outside the jet, since the mean flow velocity is zero there. Within the jet, by defining 
$\beta =\sqrt {M_{-}^{2}-1}$, we can rewrite (2.14) to be
\begin{equation} \frac{\partial^{2}p_{m-}}{\partial y^{\prime2}}-\beta^{2} \frac{\partial^{2}p_{m-}}{\partial x^{\prime2}}=0. \end{equation}
We see that the linear pressure field within the jet induced by weak shock waves satisfies the wave equation. Such an equation can admit many solutions subject to different boundary conditions; for illustration purposes, Kerschen & Cain (Reference Kerschen and Cain1995) used a step function for their single planar vortex sheet. In order to have a much more realistic shock cell structure, we use Pack's model (Pack Reference Pack1950) in this paper. Note that we use the jet height 
$D$ here, instead of the nozzle height, to non-dimensionalize the streamwise and cross-flow coordinates, i.e. 
$x=x^\prime /D$, 
$y=y^\prime /D$. The velocity potential perturbation within the jet induced by the shock waves can be written as (Pack Reference Pack1950; Tam Reference Tam1972)
\begin{equation} \phi=\sum_{j=1}A_{j}\cos(\beta a_{j} y)\sin(a_{j} x), \end{equation}
where 
$j$ represents the 
$j$th mode of the shock wave. Two parameters 
$A_{j}$ and 
$a_{j}$ in the equation above are
$$\begin{gather} A_{j}=({-}1)^{j}\frac{4\beta}{{\rm \pi}^{2}}\frac{\mathcal{U}}{(2j-1)^2}, \end{gather}$$
$$\begin{gather}a_{j}=\frac{(2j-1){\rm \pi}}{\beta}. \end{gather}$$
The constant 
$\mathcal {U}$ in the above equation is defined by 
$\mathcal {U}=U_{1}-U.$ We see from (2.17) that the amplitude of this potential function decreases quickly as the mode number increases. In light of the linearity of the model, it is convenient to consider each mode separately. In the present study, we focus on the leading-order mode. Higher-order terms can be easily included at a later stage should necessity arise. For the leading-order mode, the shock wave is periodically distributed along the streamwise direction with a shock spacing 
$s=2{\rm \pi} /a_{1}$. In what follows, the subscripts in parameters 
$A_{1}$ and 
$a_{1}$ are omitted for clarity.
The corresponding velocity, pressure and vortex-sheet deflection at the boundary of the jet flow are shown in Appendix A.
2.3. The instability waves
Similar to the derivation of the shock equations, the governing equations for the instability waves can be obtained by collecting the 
$O(\epsilon )$ terms, i.e.
$$\begin{gather} \frac{\mathrm{D}_{0}p_{v}}{\mathrm{D} t^\prime}+\rho_{0}c_{0}^{2} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{v}=0, \end{gather}$$
$$\begin{gather}\rho_{0}\frac{\mathrm{D}_{0}\boldsymbol{u}_{v}}{\mathrm{D} t^\prime}={-}\boldsymbol{\nabla} p_{v}. \end{gather}$$
Considering the 
$O(\epsilon )$ terms in the boundary conditions shown in (2.5), (2.6) and (2.7), we see that the two matching conditions can be linearized to the dynamic and kinematic conditions on 
$y=\pm D/2$, i.e.
$$\begin{gather} p_{v+}=p_{v-}, \end{gather}$$
$$\begin{gather}v_{v+}=\frac{\partial h^\prime_{v}}{\partial t^\prime},\quad v_{v-}=\frac{\partial h^\prime_{v}}{\partial t^\prime}+U\frac{\partial h^\prime_{v}}{\partial x^\prime}, \end{gather}$$
where 
$h^\prime _{v}$ denotes the disturbed height of the vortex sheets due to the instability waves. Since the initial base flow is irrotational and inviscid both inside and outside of the vortex sheet, the linear perturbation can be expressed as a velocity potential 
$\phi _{v}$. The continuity equation (2.19) and the momentum equation (2.20) can be combined to yield
\begin{equation} \boldsymbol{\nabla}^{2}\phi_{v}-\frac{1}{c_{0}^{2}}\frac{\mathrm{D}_{0}^{2}\phi_{v}}{\mathrm{D} t^{\prime2}}=0. \end{equation}Similarly the dynamic boundary condition reduces to
\begin{equation} \rho_{0+}\frac{\partial \phi_{v+}}{\partial t^\prime}=\rho_{0-} \left[\frac{\partial\phi_{v-}}{\partial t^\prime}+U\frac{\partial \phi_{v-}}{\partial x^\prime}\right]. \end{equation}For the kinematic boundary condition, the two equations shown in (2.22a,b) can be combined to yield
\begin{equation} \left(\frac{\partial}{\partial t^\prime}+U\frac{\partial}{\partial x^\prime}\right) \frac{\partial \phi_{v+}}{\partial y^\prime}=\frac{\partial^{2}\phi_{v-}}{\partial t^{'} \partial y^\prime}. \end{equation}
With temporal and spatial harmonic assumptions, the perturbations induced by the instability waves have the form of 
$\exp ({{\mathrm i}(\alpha ^\prime x^\prime -\omega ^\prime t^\prime )})$, where the spatial wavenumber 
$\alpha ^\prime$ is complex and the eigenvalue with a negative imaginary part represents instability. We use the velocity of the fully expanded jet flow 
$U$ to non-dimensionalize other variables; for instance, the non-dimensional time, frequency and wavenumber are 
$t=t^\prime U/D, \omega =\omega ^\prime D/U$ and 
$\alpha =D\alpha ^\prime$, respectively.
Combining the governing equations and boundary conditions, and noticing that the base flow outside the jet is zero, we find that the velocity potential can be expressed as
\begin{equation} \phi=U D \exp({{\rm
i}(\alpha x-\omega t)})\times \begin{cases}
\dfrac{1}{M_{+}^{2}}\,\mathrm{e}^{{-}m_{+}y}, &
y>\dfrac{1}{2} \\
\dfrac{1}{M_{-}^{2}}\dfrac{\omega}{\omega-\alpha}
(k_{1}\,\mathrm{e}^{m_{-}y}+k_{2}\,\mathrm{e}^{{-}m_{-}y}),
& y\leq \left|\dfrac{1}{2}\right|\\
\dfrac{1}{M_{+}^{2}}k_{3}\,\mathrm{e}^{m_{+}y}, &
y<{-}\dfrac{1}{2}, \end{cases}
\end{equation}
where 
$m_{+}=\sqrt {\alpha ^{2}-\omega ^{2}M_{+}^{2}}$ and 
$m_{-}=\sqrt {\alpha ^{2}-M_{-}^{2}(\omega -\alpha )^{2}}$. The branch cut is chosen such that the real part of 
$m_{+}$ is positive. Here 
$k_{1}, k_{2}, k_{3}$ are undetermined coefficients. It can be seen that both antisymmetric and symmetric modes can exist, corresponding to 
$k_{3}=- 1$ and 
$k_{3}=1$, respectively. Using the dynamic and kinematic boundary conditions, the dispersion relations can be found to be
\begin{equation} \mathrm{e}^{2m_{-}}=\frac{(\omega^{2}m_{-}/M_{-}^{2}- (\omega-\alpha)^2m_{+}/M_{+}^{2})^{2}}{(\omega^{2}m_{-}/M_{-}^{2}+ (\omega-\alpha)^2m_{+}/M_{+}^{2})^{2}}. \end{equation}For the symmetric mode, this reduces to
\begin{equation} \tanh\left(\frac{m_{-}}{2}\right) \frac{\omega^{2} m_{-}}{M_{-}^{2}(\omega-\alpha)^{2}}+\frac{m_{+}}{M_{+}^{2}}=0, \end{equation}where the parameters attain the following values:
\begin{equation} k_{1}=k_{2}=\frac{\exp\left({-\dfrac{1}{2}m_{+}}\right)}{2\cosh \left(\dfrac{1}{2}m_{-}\right)},\quad k_{3}=1. \end{equation}For the antisymmetric mode, the dispersion relationship reduces to
\begin{equation} \frac{\omega^{2} m_{-}}{M_{-}^{2}(\omega-\alpha)^{2}}+ \tanh\left(\frac{m_{-}}{2}\right)\frac{m_{+}}{M_{+}^{2}}=0, \end{equation}where the three parameters take the values of
\begin{equation} k_{1}={-}k_{2}=\frac{\exp\left({-\dfrac{1}{2}m_{+}}\right)}{2\sinh \left(\dfrac{1}{2}m_{-}\right)},\quad k_{3}={-}1. \end{equation}The corresponding pressure, velocity and deflection of the jet boundary due to the instability waves are shown in Appendix B. We can see that the deflections generated by the instability wave at the upper and lower boundaries are symmetric and antisymmetric for the symmetric and antisymmetric modes, respectively, as would be expected.
Experiments found that rectangular jets are capable of sustaining both symmetric and antisymmetric oscillation modes (Suda et al. Reference Suda, Manning and Kaji1993; Kaji & Nishijima Reference Kaji and Nishijima1996), which can be directly linked to the instability of the jet. In our analysis, both symmetric and antisymmetric instability modes can be considered. But for jet flows from high-aspect-ratio rectangular nozzles, the flapping mode is dominant (Edgington-Mitchell Reference Edgington-Mitchell2019) and the problem can be approximated by a two-dimensional theory. So in what follows only the antisymmetric mode of instability waves is considered.
2.4. The interaction between shock and instability waves
Having obtained the shock and instability waves, we are now in a position to consider the 
$O(\delta \epsilon )$ terms in the governing equations. When a perfect gas is assumed (
$\rho _0 c_0^{2}=\gamma p_0$), the continuity and momentum equations can be expressed as
$$\begin{gather} \frac{\mathrm{D}_{0}p_{i}}{\mathrm{D} t^\prime}+\rho_{0}c_{0}^{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{i}={-} \left[\boldsymbol{u}_m\boldsymbol{\cdot}\boldsymbol{\nabla} p_{v}+\boldsymbol{u}_{v}\boldsymbol{\cdot}\boldsymbol{\nabla} p_{m}+\gamma p_{m}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_v+\gamma p_{v}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{m}\right], \end{gather}$$
$$\begin{gather}\rho_{0}\frac{\mathrm{D}_{0}\boldsymbol{u}_{i}}{\mathrm{D}t^\prime}+\boldsymbol{\nabla} p_{i}={-}\left[\rho_{0}(\boldsymbol{u}_m\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_{v}+\boldsymbol{u}_{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_{m})+\rho_{m}\frac{\mathrm{D}_{0}\boldsymbol{u}_v}{{\rm D} t^\prime}+\rho_{v}\frac{\mathrm{D}_{0}\boldsymbol{u}_m}{{\rm D} t^\prime}\right]. \end{gather}$$Substituting the momentum equations for the shock wave, i.e. (2.13), and the instability wave, i.e. (2.20), into (2.33), we have
\begin{equation} \rho_{0}\frac{\mathrm{D}_{0}\boldsymbol{u}_{i}}{\mathrm{D} t^\prime}+\boldsymbol{\nabla} p_{i}={-}\left[\rho_{0}\boldsymbol{\nabla}(\boldsymbol{u}_{m}\boldsymbol{\cdot}\boldsymbol{u}_{v})+ \frac{1}{\rho_{0}c_{0}^{2}}\boldsymbol{\nabla}(p_{m}p_{v})\right]. \end{equation}
The interaction field is also irrotational, i.e. 
$u_{i}=\boldsymbol {\nabla } \phi _{i}$. Then (2.34) can be integrated to obtain
\begin{equation} p_{i}={-}\rho_{0}\frac{\mathrm{D}_{0}\phi_{i}}{\mathrm{D} t^\prime}-\rho_{0}\boldsymbol{u}_{m}\boldsymbol{\cdot}\boldsymbol{u}_{v}+ \frac{1}{\rho_{0}c_{0}^{2}}p_{m}p_{v}. \end{equation}
Combining (2.32) and (2.35), and considering the momentum and continuity equations for the 
$O(\alpha )$ and 
$O(\epsilon )$ terms, we find that 
$\phi _{i}$ satisfies the following inhomogeneous wave equation:
\begin{align} \boldsymbol{\nabla}^{2}\phi_{i}-\frac{1}{c_{0}^{2}}\frac{\mathrm{D}_{0}^{2}\phi_{i}}{\mathrm{D} t^{\prime2}}=\frac{-1}{\rho_{0}c_{0}^{2}}[2(\boldsymbol{u}_m\boldsymbol{\cdot}\boldsymbol{\nabla} p_{v}+\boldsymbol{u}_v\boldsymbol{\cdot}\boldsymbol{\nabla} p_{m})+(\gamma-1)(p_{m}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{v}+p_{v} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{m})]. \end{align}
Next we consider the continuity of the pressure across the vortex sheet. Similar to Kerschen & Cain (Reference Kerschen and Cain1995), on the two boundaries 
$y^\prime =1/2D$ and 
$y^\prime =-1/2D$, the dynamic boundary condition reduces to
\begin{equation} p_{i+}+h_{m}^\prime\frac{\partial p_{v+}}{\partial y^\prime}= p_{i-}+h_{m}^\prime \frac{\partial p_{v-}}{\partial y^\prime}+h_{v}^\prime\frac{\partial p_{m-}}{\partial y^\prime}, \end{equation}while the kinematic boundary condition requires
$$\begin{gather} v_{i+}+\frac{\partial v_{v+}}{\partial y^\prime}h^\prime_{m}= \frac{\partial h_{i}^\prime}{\partial t^\prime}+u_{v+}\frac{\partial h_{m}^\prime}{\partial x^\prime}, \end{gather}$$
$$\begin{gather}v_{i-}+\frac{\partial v_{v-}}{\partial y^\prime}h^\prime_{m}+\frac{\partial v_{m-}}{\partial y^\prime}h^\prime_{v}=\frac{\partial h_{i}^\prime}{\partial t^\prime}+U\frac{\partial h^\prime_{i}}{\partial x^\prime}+u_{v-}\frac{\partial h_{m}^\prime}{\partial x^\prime}+u_{m-}\frac{\partial h_{v}^\prime}{\partial x^\prime}. \end{gather}$$
Note that outside the jet flow, there is no perturbation due to shock waves, therefore the right-hand side of (2.36) vanishes. This equation degenerates to a homogeneous wave equation. Using the same 
$U$ and 
$D$ to non-dimensionalize the velocity potential, we obtain
$$\begin{gather} \phi_{i+}=2 U D g_{+}\exp({-\mathrm{i}\omega t}), \end{gather}$$
$$\begin{gather}\phi_{i-}=2 U D g_{-}\exp({-\mathrm{i}\omega t}), \end{gather}$$
where 
$g_\pm$ are the non-dimensionalized potential functions. Outside the jet flow, the base flow is uniformly zero, and (2.36) reduces to
\begin{equation} \frac{\partial^{2}g_{+}}{\partial x^{2}}+ \frac{\partial^{2}g_{+}}{\partial y^{2}}+\omega^{2}M_{+}^{2}g_{+}=0. \end{equation}
Inside the jet flow, 
$g_{-}$ satisfies
\begin{align} \frac{\partial^{2}g_{-}}{\partial x^{2}}+\frac{\partial^{2}g_{-}}{\partial y^{2}}+M_{-}^{2}(\omega+\mathrm{i}\frac{\partial}{\partial x})^{2}g_{-} &= 2k_{1}\,\mathrm{e}^{\mathrm{i} \alpha x} \frac{A}{U}[\sinh(\zeta_{1}y)(B_{1}\cos(a x)+B_{2}\sin(a x))\nonumber\\ &\quad + \sinh(\zeta_{2}y)(B_{3}\cos(a x)+B_{4}\sin(a x)], \end{align}where
\begin{equation} \left.\begin{gathered} \zeta_{1}=m_{-}+\mathrm{i} a\beta,\\ \zeta_{2}=m_{-}- \mathrm{i}a\beta, \end{gathered}\right\} \end{equation}
\begin{equation} \left.\begin{gathered} B_{1}=\frac{a\omega}{2}\left(\alpha+\mathrm{i}\frac{a\beta m_{-}}{\omega-\alpha}- \frac{\gamma-1}{2}(\omega-\alpha)M_{-}^{2}\right),\\ B_{2}=\frac{a\omega}{2}\left(m_{-}\beta-\mathrm{i}\frac{\alpha a}{\omega-\alpha}+ \frac{\gamma-1}{2}\mathrm{i}a M_{-}^{2} \right),\\ B_{3}=\frac{a\omega}{2}\left(\alpha-\mathrm{i}\frac{a\beta m_{-}}{\omega-\alpha}- \frac{\gamma-1}{2}(\omega-\alpha)M_{-}^{2}\right),\\ B_{4}=\frac{a\omega}{2}\left({-}m_{-}\beta-\mathrm{i}\frac{\alpha a}{\omega-\alpha}+\frac{\gamma-1}{2}\mathrm{i}a M_{-}^{2}\right). \end{gathered}\right\} \end{equation}Besides, the two boundary conditions can be reorganized as
\begin{align} \omega g_{+}-\frac{M_{-}^{2}}{M_{+}^{2}}(\omega+\mathrm{i}\frac{\partial}{\partial x})g_{-}&={\pm}\frac{A\beta}{2UM_{+}^{2}}\sin\left(\frac{1}{2}a\beta \right) \cos(ax)\left[\omega\left(2k_{1}m_{-}\cosh\left(\frac{1}{2}m_{-}\right) \vphantom{\left(\frac{1}{2}m_{-}\right)}\right. \vphantom{\frac{\exp\left({-\dfrac{1}{2}m_{+}}\right)}{\omega}}\right. \nonumber\\ &\quad\pm\left.\left.\vphantom{\left(\frac{1}{2}m_{-}\right)} m_{+}\exp\left({\mp\frac{1}{2}m_{+}}\right)\right) -\frac{\exp\left({-\dfrac{1}{2}m_{+}}\right)}{\omega}m_{+} \frac{M_{-}^{2}}{M_{+}^{2}}a^{2}\right]\mathrm{e}^{\mathrm{i}\alpha x} \end{align}and
\begin{align} \left(1+\frac{\mathrm{i}}{\omega}\frac{\partial}{\partial x}\right)\frac{\partial g_{+}}{\partial y}-\frac{\partial g_{-}}{\partial y}&=\frac{A}{2 U} \left[{\mathrm{i}}a\beta\sin(ax)\left [b_{1}\sin\left({\pm}\frac{1}{2}a\beta\right)+c\beta \cos\left(\frac{1}{2}a\beta\right)\right]\right. \nonumber\\ &\quad +\left.\cos(a x)\left[b_{2}\sin\left({\pm}\frac{1}{2}a\beta\right)-c \alpha\cos\left(\frac{1}{2}a\beta\right)\right]\right]\mathrm{e}^{\mathrm{i} \alpha x}, \end{align}
where the upper and lower of signs of 
$\pm$ and 
$\mp$ correspond to the matching conditions on 
$y=1/2$ and 
$y=-1/2$, respectively, and
$$\begin{gather} b_{1}={\pm}\frac{2\alpha k_{1}\sinh\left(\dfrac{1}{2}m_{-}\right)}{M_{+}^{2}\omega} \left[(\alpha-\omega)+\frac{m_{+}^{2}}{\alpha}\right]\pm \frac{\alpha}{M_{-}^{2}}\frac{\omega}{\omega-\alpha}\exp\left({\mp\frac{1}{2}m_{+}}\right), \end{gather}$$
$$\begin{gather}b_{2}={\pm}\frac{2\beta k_{1}\sinh\left(\dfrac{1}{2}m_{-}\right)}{M_{+}^{2}\omega}[(\alpha-\omega)m_{+}^{2}+\alpha a^{2}]\pm m_{-}^{2}\frac{\beta}{M_{-}^{2}}\frac{\omega}{\omega-\alpha}\exp\left({\mp\frac{1}{2}m_{+}}\right), \end{gather}$$
$$\begin{gather}c=\frac{a}{\omega}\frac{m_{+}}{M_{+}^{2}}\exp\left({-\frac{1}{2}m_{+}}\right). \end{gather}$$
To obtain the solution to these equations, Fourier transform is used. The Fourier transforms 
$G_{\pm }(\lambda,y)$ are defined as
\begin{equation} G_{{\pm}}(\lambda,y)=\int_{-\infty}^{+\infty} g_{{\pm}}(x,y)\,\mathrm{e}^{\mathrm{i} \lambda x}\,\mathrm{d}\kern0.7pt x, \end{equation}
where 
$\lambda$ is the wavenumber in the streamwise direction. Outside the jet flow, it is easy to find that 
$G$ satisfies
\begin{equation} G_{+}(\lambda,y)=\begin{cases} D_{1}(\lambda)\,\mathrm{e}^{-\gamma_{+}y}, & y>1/2,\\ D_{4}(\lambda)\,\mathrm{e}^{\gamma_{+}y}, & y<{-}1/2, \end{cases} \end{equation}where
\begin{equation} \gamma_{+}(\lambda)=\sqrt{\lambda^{2}-\omega^{2}M_{+}^{2}} \end{equation}
and 
$D_{1}$ and 
$D_{4}$ are two undetermined coefficients related to 
$\lambda$.
When Fourier transform is used to the source term of (2.43), we can simplify the problem by noting the periodicity of the functions 
$\cos (ax)$ and 
$\sin (a x)$. For example,
\begin{align} & \int_{0}^{\infty} \sin(a x)\exp({\mathrm{i}(\lambda+\alpha) x})\,\mathrm{d}\kern0.7pt x= \sum_{j=0}^\infty \int_{2j{\rm \pi}/a}^{2(j+1){\rm \pi}/a} \sin(a x)\exp({\mathrm{i}(\lambda+\alpha) x})\,\mathrm{d}\kern0.7pt x \nonumber\\ &\quad = \left(\sum_{j=0}^\infty \exp\left({\mathrm{i}\left(\alpha+\lambda\right) 2 j {\rm \pi}/ a}\right)\right) \frac{1}{2}\left(\frac{1}{\alpha+\lambda+a}-\frac{1}{\alpha+\lambda-a}\right) \nonumber\\ &\qquad \times\left(-\exp\left({\mathrm{i}(\alpha+\lambda)\frac{2{\rm \pi}}{a}}\right)+1\right). \end{align}
Clearly, a convergence problem arises when 
$j \to +\infty$ since 
$\alpha$ has a negative imaginary part. However, this is a difficulty resulting from linearization, not from inherent difficulties in the flow physics and its modelling. As mentioned earlier, the interaction between shock and instability waves occurs when the instability waves grow to be of sufficient amplitude, at which place nonlinear/linear saturation or even decay begins to take place. Considering that both the instability and shock waves start to decay further downstream of the jet (Cohen & Wygnanski Reference Cohen and Wygnanski1987; Jordan & Colonius Reference Jordan and Colonius2013), the effective interaction only takes place within a limited interval spanning several shock cells (Suda et al. Reference Suda, Manning and Kaji1993; Gao & Li Reference Gao and Li2010; Malla & Gutmark Reference Malla and Gutmark2017). In other words, the summation only involves a finite number of terms and therefore the convergence problem does not occur in realistic jets. In light of this, it is reasonable to only focus on a limited number of shock cells in this paper, e.g. from the third to the fifth according to previous studies (Norum Reference Norum1983; Panda Reference Panda1999; Tam & Parrish Reference Tam and Parrish2014; Mercier, Castelain & Bailly Reference Mercier, Castelain and Bailly2017). Furthermore, note that the right-hand sides of (2.46), (2.47) and (2.43) all have the common factor 
$(\sum \exp ({\mathrm {i}( \alpha +\lambda ) 2 j {\rm \pi}/ a}))$ after the Fourier transform. Hence they can be collected and accounted for later due to linearity of the equation and the two boundary conditions. This means that we only need to consider the interaction within one shock cell, as shown in figure 2, and the total interaction field would be a simple linear combination from a number of shock cells. In this way, the effective integration interval when Fourier transform is applied to both the governing equation and boundary conditions is limited to be within one shock cell. Physically, this corresponds to the effective sound generated by the interaction between the instability waves and one single shock cell, similar to Powell's idea of treating the effective sound as that of monopoles.
Figure 2. Schematic of an effective source term located within one single shock cell. The total sound field is equivalent to a linear superimposition of the results from a number of shock cells.
Note although we limit the integration interval to be within one shock cell, we do not imply that such an effective source is physically localized. The source term on the right-hand side of (2.36) has a periodic nature by construction, but the bounds of integration may be across one or several shock cells due to the linearity of (2.36). This is equivalent to decomposing the problem into several subproblems, each of which has an effective noise source within one shock cell. The overall sound is a linear combination of the solutions to these subproblems. In the rest of this paper, we focus on examining the characteristics of such an effective sound source first and then discuss and compare the total sound from a number of these sources with experiments and numerical simulations.
Let us define
\begin{align} \mathcal{I}_s(\lambda) &= \int_{0}^{{2{\rm \pi}}/{a}}\sin(ax)\exp({\mathrm{i}(\lambda+\alpha)x})\,\mathrm{d}\kern0.7pt x \nonumber\\ &=\frac{1}{2}\left(\frac{1}{\alpha+\lambda+a}-\frac{1}{\alpha+\lambda-a}\right) \left(-\exp\left({\mathrm{i}(\alpha+\lambda)\frac{2{\rm \pi}}{a}}\right)+1\right), \end{align}
\begin{align} \mathcal{I}_c(\lambda) &= \int_{0}^{{2{\rm \pi}}/{a}}\cos(ax)\exp({\mathrm{i}(\lambda+\alpha)x})\,\mathrm{d}\kern0.7pt x \nonumber\\ &=\frac{\mathrm{i}}{2}\left(\frac{1}{\alpha+\lambda+a}+ \frac{1}{\alpha+\lambda-a}\right)\left(-\exp\left({\mathrm{i}(\alpha+\lambda)\frac{2{\rm \pi}}{a}}\right)+1\right), \end{align}with which the inhomogeneous equation can be written as
\begin{equation} \frac{\partial^{2}G_{-}}{\partial y^{2}}-\gamma_{-}^{2}G_{-} = 2k_{1}\frac{A}{U}\left[\sinh(\zeta_{1}y)(B_{1}\mathcal{I}_{c}+B_{2}\mathcal{I}_{s})+ \sinh(\zeta_{2}y)(B_{3}\mathcal{I}_{c}+B_{4}\mathcal{I}_{s})\right], \end{equation}where
\begin{equation} \gamma_{-}=\sqrt{\lambda^{2}-M_{-}^{2}(\omega+\lambda)^{2}}. \end{equation}Equation (2.58) is equivalent to
\begin{equation} \gamma_{-}={-}\mathrm{i}\beta\sqrt{(\lambda-M_{1})(\lambda-M_{2})}, \end{equation}where
\begin{equation} M_{1}=\frac{-M_{-}\omega}{M_{-}+1},\quad M_{2}=\frac{-M_{-}\omega}{M_{-}-1}. \end{equation}
The branch cuts passing 
$\lambda =M_{1}$ and 
$\lambda =M_{2}$ extend to the lower half-plane, as illustrated in figure 3. The function 
$G_{-}(\lambda,y)$ can be divided into two parts, a particular solution, 
$G^{p}(\lambda,y)$, and a complementary solution, 
$G^{c}(\lambda,y)$, i.e.
\begin{equation} G_{-}(\lambda,y)=G^{p}(\lambda,y)+G^{c}(\lambda,y). \end{equation}The particular solution can be calculated analytically, i.e.
\begin{equation} G^{p}(\lambda,y) =2k_{1}\frac{A}{U}\left[\frac{\sinh(\zeta_{1}y)}{\zeta_{1}^{2}-\gamma_{-}^2} (B_{1}\mathcal{I}_{c}+B_{2}\mathcal{I}_{s})+ \frac{\sinh(\zeta_{2}y)}{\zeta_{2}^{2}-\gamma_{-}^2} (B_{3}\mathcal{I}_{c}+B_{4}\mathcal{I}_{s})\right], \end{equation}and the complementary solution can be found to be
\begin{equation} G^{c}(\lambda,y)=D_{2}(\lambda)\,\mathrm{e}^{\gamma_{-}y}+D_{3}(\lambda)\,\mathrm{e}^{-\gamma_{-}y}, \end{equation}
where 
$D_{2}$ and 
$D_{3}$ are two undetermined coefficients.
Figure 3. The branch points, branch cuts and integral path in complex 
$\lambda$ plane. Paths 
${P}_1$ and 
${P}_2$ are the deformed integral paths when the saddle point approaches the branch points 
$-\omega M_{+}$ and 
$\omega M_{+}$, respectively.
Applying the Fourier transform to the two boundary conditions, we can solve the undetermined coefficients. For the antisymmetric mode, we obtain
\begin{align} D_{1}(\lambda) &= \frac{A\exp\left({\dfrac{1}{2}\gamma_{+}}\right)}{2U\eta M_{+}^{2}} \left[\gamma_{-}\coth\left(\frac{\gamma_{-}}{2}\right)\mathcal{I}_{c}\sin \left({\pm}\frac{1}{2}a\beta\right) \vphantom{\frac{M_{-}^{2}}{M_{+}^{2}}}\right.\nonumber\\ &\quad \times \beta\left[\omega\left(2k_{1}m_{-}\cosh \left(\frac{1}{2}m_{-}\right)\pm m_{+}\exp\left({\mp\frac{1}{2}m_{+}}\right)\right) \vphantom{\frac{M_{-}^{2}}{M_{+}^{2}}}\right. \nonumber\\ &\quad -\left.\frac{1}{\omega}m_{+}\frac{M_{-}^{2}}{M_{+}^{2}}a^{2} \exp\left({-\frac{1}{2}m_{+}}\right)\right]+M_{-}^{2}(\omega+\lambda) \left(\gamma_{-}\coth\left(\frac{\gamma_{-}}{2}\right)\frac{2U}{A}G^{p} \left({\pm}\frac{1}{2}\right)\right. \nonumber\\ &\quad - \frac{2U}{A}{G^{p}}^{\prime}\left({\pm}\frac{1}{2}\right)-{\rm i} a \beta \mathcal{I}_{s}\left[b_{1}\sin\left({\pm}\frac{1}{2}a\beta\right)+c\beta\cos \left(\frac{1}{2}a\beta \right)\right]\nonumber\\ &\quad\left.\vphantom{\frac{M_{-}^{2}}{M_{+}^{2}}} \left. -\mathcal{I}_{c}\left[b_{2}\sin \left({\pm}\frac{1}{2}a\beta\right)-c\alpha\cos \left(\frac{1}{2}a\beta \right)\right]\right)\right], \end{align}
\begin{align} D_{2}(\lambda) &={-} \frac{A}{4\eta U}\left [\frac{\omega+\lambda}{\omega}\mathcal{I}_{c}\gamma_{+}\beta\sin \left ({\pm}\frac{1}{2}a\beta \right )\frac{1}{M_{+}^{2}} \vphantom{\frac{M_{-}^{2}}{M_{+}^{2}}}\right.\nonumber\\ &\quad \times \left[\omega\left(2k_{1}m_{-}\cosh \left (\frac{1}{2}m_{-}\right )\pm m_{+}\exp\left({\mp\frac{1}{2}m_{+}}\right)\right) \vphantom{\frac{M_{-}^{2}}{M_{+}^{2}}}\right. \nonumber\\ &\quad -\left.\frac{1}{\omega}m_{+}\frac{M_{-}^{2}}{M_{+}^{2}}a^{2} \exp\left({-\frac{1}{2}m_{+}}\right)\right]+\frac{M_{-}^{2}}{M_{+}^{2}} \frac{(\omega+\lambda)^{2}}{\omega}\gamma_{+}\frac{2U}{A}G^{p}\left({\pm}\frac{1}{2}\right) \nonumber\\ &\quad +\frac{2U}{A}{G^{p}}^{\prime}\left({\pm}\frac{1}{2}\right)+ \mathrm{i}\omega a \beta \mathcal{I}_{s}\left[b_{1}\sin \left({\pm}\frac{1}{2}a\beta\right)+c\beta\cos \left(\frac{1}{2}a\beta\right)\right] \nonumber\\ &\quad \left.\vphantom{\frac{M_{-}^{2}}{M_{+}^{2}}} +\omega\mathcal{I}_{c}\left[b_{2}\sin\left({\pm}\frac{1}{2}a\beta\right)-c\alpha \cos\left(\frac{1}{2}a\beta\right)\right]\right], \end{align}
$D_{4}(\lambda )=-D_{1}(\lambda )$ and 
$D_{3}(\lambda )=-D_{2}(\lambda )$, where 
$\eta (\lambda )$ is
\begin{equation} \eta(\lambda)=\omega{\rm coth}\left(\frac{1}{2}\gamma_{-}\right)\gamma_{-}+\frac{1}{\omega}(\omega+\lambda)^{2}\gamma_{+} \frac{M_{-}^{2}}{M_{+}^{2}}. \end{equation}
It is straightforward to verify that 
$\lambda =-\alpha$ is a simple zero for 
$\eta (\lambda )$. In fact, 
$\eta (-\alpha )=0$ corresponds to the dispersion relation (2.30). Besides, 
$\lambda =-\alpha$ is a simple zero for 
$\mathcal {I}_s(\lambda )$ and a second-order zero for 
$\mathcal {I}_c(\lambda )$, so 
$\lambda =-\alpha$ is not a pole for 
$D_{1}(\lambda )$ and 
$D_{2}(\lambda )$. It is found that 
$D_{1}(\lambda )$ and 
$D_{2}(\lambda )$ have four poles, which are
$$\begin{gather} Z_{1,2}=\frac{\omega M_{-}^2\pm\sqrt{M_{-}^{2}\omega^{2}+(1-M_{-}^{2})\zeta_{1}^{2}}}{1-M_{-}^2}, \end{gather}$$
$$\begin{gather}Z_{3,4}=\frac{\omega M_{-}^2\pm\sqrt{M_{-}^{2}\omega^{2}+(1-M_{-}^{2})\zeta_{2}^{2}}}{1-M_{-}^2}. \end{gather}$$The following inverse Fourier transform:
\begin{equation} g_{+}(x,y)=\frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty}D_{1}(\lambda)\exp({-(\mathrm{i}\lambda x+\gamma_{+}y)})\,\mathrm{d}\lambda \end{equation}
yields 
$g_{+}$. The integration path is near the real axis of 
$\lambda$, as illustrated in figure 3. The integration path is indented to pass above the poles at 
$\lambda =Z_{1,4}$, as illustrated in figure 3, in accordance with the causality argument (Briggs Reference Briggs1964). Because the real part of 
$\gamma _{+}$ should be positive when 
$|\lambda |\rightarrow \infty$ along the integration path, the branch cuts of 
$\gamma _{+}$ passing the branch points 
$\lambda =\pm \omega M_{+}$ are chosen to extend to the upper and lower half-plane, respectively, as shown in figure 3. The branch points of 
$\gamma _{-}$, i.e. 
$M_{1}$ and 
$M_{2}$, are on the negative real 
$\lambda$ axis. The branch cuts are chosen to extend down to the lower half-plane so as not to cross the integration path. Using the steepest descent method, and noting that the saddle point is located at 
$\lambda =-M_{+}\omega \cos \theta$, where 
$\theta =\arctan (y/x)$ represents the observer angle, we can express 
$g_{+}$ as a function of radial distance 
$r$ and 
$\theta$ in the far field (
$r\gg 1$), i.e.
\begin{equation} g_{+}(r,\theta)=\frac{\sqrt{M_{+}\omega}}{\sqrt{2{\rm \pi}}}D_{1}({-}M_{+}\omega\cos \theta)\sin \theta \frac{\exp({\mathrm{i}\omega(M_{+}r-{\rm \pi}/4)})}{\sqrt{r}}+O(r^{{-}3/2}), \end{equation}
and with (2.35) and (2.40), the corresponding pressure perturbation (non-dimensionalized by 
$\sqrt {2/{\rm \pi} }\rho _{0+}U^2$) can be expressed as
\begin{align} p_{+}(r,\theta)=\mathrm{i}\sqrt{M_{+}}\omega^{{3}/{2}}D_{1}({-}M_{+}\omega\cos \theta){\sin} \theta \frac{\exp({\mathrm{i}\omega(M_{+}r-t-{\rm \pi}/4)})}{\sqrt{r}}+O(r^{{-}3/2}). \end{align}
Note that the saddle point moves between 
$-\omega M_{+}$ and 
$\omega M_{+}$ as 
$\theta$ changes from 0 to 
${\rm \pi}$, and the integral path is forbidden to pass through the branch cut. So when 
$\theta ={\rm \pi}$ and 
$0$, the corresponding integral path is deformed along the branch cut and wraps the branch point as shown by 
${P}_1$ and 
${P}_2$ in figure 3, respectively. It is similar when the steepest descent path passes through the branch cut at 
$M_{1}$ and 
$M_{2}$, in which case the integral path needs to be adjusted to avoid the branch cut. Note that when the poles cross the steepest descent path, care must be taken regarding the residue contribution.
As can be seen from (2.70), when the observer angle 
$\theta ={\rm \pi}$, the potential function 
$g_{+}(r,\theta )$ reduces to a high-order term 
$O(r^{-3/2})$ if 
$|D_{1}(-\omega M_{+}\cos \theta )|$ is bounded. This implies that sound waves propagating in this direction decay rapidly and nearly vanish in the far field as 
$r\rightarrow \infty$. This leads to an important feature of the sound directivity that will become clear in the rest of this paper.
3. Results and discussion
The sound field due to the interaction between shock and instability waves is shown in this section. In the linear stability analysis, the spatial wavenumber 
$\alpha$ is the central parameter determining the characteristics of the instability waves. The dispersion relation calculated from (2.30) is shown in figure 4. The antisymmetric mode is considered. We see that both the wavenumber and growth rates increase as 
$\omega$ increases. These are well-established results in the linear stability analysis. Due to a negative imaginary part, the instability waves grow exponentially downstream of the jet flow, and subsequently interact with shock cell structures.
Figure 4. The real and imaginary parts of the solution of the spatial wavenumber 
$\alpha$ calculated from the dispersion relation (2.30). The antisymmetric mode is considered. Note that the imaginary part of the wavenumber 
$\alpha$ is negative, and its absolute value is plotted instead.
In what follows, we first examine the screech frequency prediction using the present model. Sound propagating at the observer angle 
$\theta =150^\circ$ is used to verify the far-field approximation before the directivity patterns of the fundamental tone and its harmonics are shown. Finally, we examine the near-field pressure fluctuation and discuss the noise generation mechanism.
3.1. The screech frequency
Powell (Reference Powell1953b) proposed a model to predict the screech frequency by assuming a constructive interference in the upstream direction 
$\theta =180^\circ$. Following Powell's idea, the screech frequency and its harmonics can be calculated using the present model. Note that similar frequency predictions have been investigated in earlier studies; nevertheless, it is included here as a validation of the model. The shock cell spacing 
$s$ satisfies
\begin{equation} s/D=\frac{2{\rm \pi}}{a}=2\sqrt{M_{-}^{2}-1}. \end{equation}Tam (Reference Tam1986) showed that in the jet flow from a rectangular nozzle, the shock cell spacing satisfies
\begin{equation} s/D=\frac{2\sqrt{M_{-}^{2}-1}}{\sqrt{1+\left(\dfrac{D}{b}\right)^{2}}}, \end{equation}
where 
$b$ is the width of the fully expanded jet flow. When the rectangular jet is of high aspect ratio 
$(D/b\ll 1)$, we can see that (3.2) reduces to (3.1). Note that (3.1) may not be able to predict the shock spacing accurately when the magnitude of overexpansion/underexpansion increases (Gao & Li Reference Gao and Li2009), while a correct representation of shock structures plays a dominant role in predicting the screech frequency, especially in circular jets (Nogueira et al. Reference Nogueira, Jaunet, Mancinelli, Jordan and Edgington-Mitchell2022a). Nevertheless, in this planar model, the jet is assumed to be slightly imperfectly expanded, and previous studies (Tam Reference Tam1986) showed good agreement between the prediction and the experimental data. Therefore, we choose (3.1) to predict the shock spacing.
The convection velocity of instability waves is widely believed to be proportional to the velocity of the fully expanded jet flow, i.e.
\begin{equation} U_{c}=\kappa U_{-}, \end{equation}
where 
$\kappa$ is usually taken to be 
$0.7$. Equations (1.3), (3.1) and (3.3) can be combined to yield the non-dimensionalized angular frequency of the sound wave:
\begin{equation} \omega=\frac{a m \kappa }{M_c +1}, \end{equation}
where 
$m=1, 2, 3, 4,\ldots$ correspond to the fundamental frequency and the first, second and higher harmonics, respectively. Considering that 
$M_c=U_c/a_\infty$, where 
$a_\infty$ is the speed of sound in the free stream, for a cold jet, 
$M_c$ can be calculated by
\begin{equation} M_c=\frac{\kappa M_{-}}{\sqrt{1+\dfrac{\gamma-1}{2}M_{-}^2}}. \end{equation}
With (3.5), 
$\omega$ can be readily calculated via (3.4). This formula is consistent with that derived by Tam (Reference Tam1986). A comparison between the measured fundamental screech frequency by Powell (Reference Powell1953b) and that predicted by (3.4) is shown in figure 5. As can be seen, good agreement is achieved. Equations (3.4) and (3.5) show that when the jet operating condition is known, the screech frequency can be readily calculated. The operating conditions and frequencies calculated in this way, as shown in table 1, are used in subsequent sections.
Figure 5. Comparison of the fundamental screech frequency between Powell's experiment (Powell Reference Powell1953b) and the model's prediction.
Table 1. The operating conditions and calculated frequencies of the screech tones, where 
$\omega _{1}$ denotes the angular frequency of the fundamental tone, 
$\omega _{2}$ the first harmonic and 
$\omega _{3}$ the second harmonic.
Note that this paper does not attempt to model the entire feedback loop, and the reason we include a frequency prediction is mainly to validate the model. Thus, although the feedback theory proposed by Powell is used, it is only used to predict the frequency using (3.4). Our focus in this paper is to model the interaction between the shock and instability waves. It is worth noting that a new feedback mechanism for circular jets has been proposed in a number of recent papers (Edgington-Mitchell et al. Reference Edgington-Mitchell, Jaunet, Jordan, Towne, Soria and Honnery2018; Gojon et al. Reference Gojon, Bogey and Mihaescu2018; Li et al. Reference Li, Zhang, Hao and He2020; Mancinelli et al. Reference Mancinelli, Jaunet, Jordan and Towne2021). The present paper, however, focuses on a two-dimensional jet, and it is not yet clear what role the guided jet modes play in this case. Besides, even for circular jets, it is known that the convection velocity of the upstream-travelling jet guided mode is very close to the speed of sound. Therefore, it can be expected that little change in the frequency prediction would occur even if the guided jet mode is taken as the closure mechanism.
3.2. Directivity of the sound field
The distinct directivity pattern of jet screech is perhaps one of its most important features and has been well reported in various experiments. In this section, we aim to predict the noise directivity using the model developed in § 2. As mentioned in § 1, this paper concerns the acoustic emission due to the interaction between shock and instability waves, and therefore does not consider the entire feedback process of screech. However, the existence of the feedback would not alter the fact that the noise is generated due to the shock–instability interaction. Therefore, provided the screech frequency is specified, the present model can be used to compare with the screech directivity.
Before being used to study the directivity of the resulting sound, (2.71) is verified by numerically integrating (2.69) at 
$\theta =150^\circ$. The comparison between the prediction using (2.71) and (2.69) is shown in figure 6. The sound pressure level (SPL) is defined to be
\begin{equation} \mathrm{SPL}=10\log_{10}\frac{|p_{+}|^2}{|p_{{ref}}|^2}, \end{equation}
where the reference pressure 
$p_{ref}=2\times 10^{-5}$. We see that good agreement is achieved when 
$r$ is beyond 
$5$, where the difference between the two methods is within 1 dB. When 
$r\geq 20$, the difference reduces to 0.2 dB. Consequently 
$r=5$, in which case 
$kr\approx 7.9$, may be used to approximately separate the acoustic near and far field (Arndt, Long & Glauser Reference Arndt, Long and Glauser1997).
Figure 6. Comparison of the SPL from the far-field approximation and that from numerical integration (
$\theta =150^\circ$). The Mach number of the fully expanded jet is 1.5. The origin 
$r=0$ represents the beginning of the first shock cell, and 
$r$ is the non-dimensionalized radial distance.
Whether the sound source of jet screech is spatially localized or distributed along the jet flow is still open to debate. As mentioned above, Powell proposed the monopole array theory to predict the directivity patterns. Many other researchers (Kaji & Nishijima Reference Kaji and Nishijima1996; Raman Reference Raman1997; Malla & Gutmark Reference Malla and Gutmark2017) also found that the fundamental screech tone was emitted from several shock cells downstream of the jet flow. However, other researchers (Walker & Thomas Reference Walker and Thomas1997; Mercier et al. Reference Mercier, Castelain and Bailly2017) observed that the screech was produced from a particular shock cell, e.g. the third or fourth one. In addition, ‘shock clapping’ (Suda et al. Reference Suda, Manning and Kaji1993) and ‘shock leakage’ (Suzuki & Lele Reference Suzuki and Lele2003) were observed at the third and fourth shock cells downstream of the jet, which were suggested to be the source of screech, particularly for higher harmonics (Semlitsch et al. Reference Semlitsch, Malla, Gutmark and Mihaescu2020). This suggests that the number of shock cells that need to be included is not clear. However, it would be interesting to compare and contrast the noise directivity patterns due to the interaction between the instability waves and various numbers of shock cells. Therefore, in what follows we examine the directivity patterns due to the interaction between the instability wave and a single shock cell, and then discuss that from several shock cells.
It is known that in realistic jets the instability waves exhibit a characteristic structure of wave packets (Lele Reference Lele2005; Jordan & Colonius Reference Jordan and Colonius2013; Wong et al. Reference Wong, Jordan and Honnery2019). The amplitude of the instability waves varies slowly within one wavelength (Suzuki & Colonius Reference Suzuki and Colonius2006), while the whole wave packet shows a Gaussian envelope (Freund Reference Freund2011), or more precisely an exponentially modified Gaussian envelope as demonstrated by recent work (Maia et al. Reference Maia, Jordan, Cavalieri and Jaunet2019). We see that the local growth rate within a wavelength is varying, but the effects of the local growth rate on the sound characteristics are not clear across the wave packet. In this paper, we also examine the effects of the local growth rate by showing results with various values of 
$\alpha _{{i}}$.
3.2.1. Directivity pattern of sound due to one-cell interaction
When the fully expanded Mach number is given, the directivity patterns of the fundamental tone and its harmonics can be calculated from (2.71). The operating condition is shown in table 1. Note that the wavenumber is obtained from (2.30), as we consider the antisymmetric mode. The directivity patterns of the fundamental tone and its first two harmonics under three different operating conditions are shown in figure 7. The SPL is defined by (3.6). Labels (1), (2) and (3) represent the results for the fully expanded jet Mach number of 1.5, 1.3 and 1.2, respectively. Figures 7(a)–7(c) show the results of the fundamental tone and the first and second harmonics, respectively. From figure 7, it is clear that the effective directivity of the fundamental tone due to a single shock cell is not that of a monopole. Instead, it consists of two lobes. One primary lobe radiates upstream, while the other radiates downstream with a weaker intensity. Although this represents the effective directivity due to one shock cell interaction, we see that it possesses some inherent directivity that resembles the total sound field measured in experiments. For example, the effective directivity of the fundamental tones accords with both numerical (Berland et al. Reference Berland, Bogey and Bailly2007) and experimental (Walker & Thomas Reference Walker and Thomas1997; Malla & Gutmark Reference Malla and Gutmark2017) results for rectangular jets of high aspect ratios, where both an upstream lobe and a downstream lobe of similar intensities appear. However, a precise match of every lobe position may not be achieved; this is expected because the prediction is only from one shock cell interaction, whereas the numerical and experimental results are from the entire screeching jet. Moreover, in all three cases, the fundamental tone is reinforced at the upstream direction, but drops quickly as 
$\theta$ approaches 
$180^\circ$. In the original paper of Kerschen & Cain (Reference Kerschen and Cain1995), the directivity pattern reaches its maximum at 
$48^\circ$ with only one lobe in the downstream direction. Little sound radiates in the upstream direction. Our present model shows that when a two-dimensional vortex sheet and a more realistic shock structure are considered, the screech directivity from a single shock cell has a major radiation lobe in the upstream direction. This result shows a better qualitative agreement with experiments. It is also interesting to note that the maximal radiation angles shown in figure 7 appear to depend on the fully expanded jet Mach number. A similar tendency was also reported in earlier experiments for round jets (Powell, Umeda & Ishii Reference Powell, Umeda and Ishii1992).
Figure 7. The directivity of sound in the far field obtained by (2.71). Here 
$r$ is fixed to be 1. Labels (1), (2) and (3) represent the results for a fully expanded jet Mach number of 1.5, 1.3 and 1.2, respectively. Results of (a) the fundamental tone, (b) the first harmonic and (c) the second harmonic. The antisymmetric mode of instability waves is taken, and the imaginary part of wavenumber 
$\alpha _{{i}}\neq 0$. In addition, 
$\mathcal {U}$ in (2.17) is taken to be 1.
Another important result is that there is a large lobe perpendicular to the jet flow for the first harmonic in all three conditions, which again resembles the total noise directivity measured in experiments (Powell Reference Powell1953b; Walker & Thomas Reference Walker and Thomas1997; Semlitsch et al. Reference Semlitsch, Malla, Gutmark and Mihaescu2020). Note that the peak radiation angle was reported to be not exactly at 
$90^\circ$, but slightly towards the downstream direction (Berland et al. Reference Berland, Bogey and Bailly2007; Kandula Reference Kandula2008; Semlitsch et al. Reference Semlitsch, Malla, Gutmark and Mihaescu2020), which is similar to the prediction of the present model. For the second harmonic, it was experimentally observed that the effective directivity pattern showed two lobes, one directed slightly upstream and one downstream (Malla & Gutmark Reference Malla and Gutmark2017); this feature is consistent with the prediction. In addition, we find that in a recent experiment conducted by Semlitsch et al. (Reference Semlitsch, Malla, Gutmark and Mihaescu2020), the main radiation angle for the second harmonic was between 
$40^\circ$ and 
$110^\circ$, while little radiation appeared around 
$90^\circ$. We see from figure 7 that this is also reflected in the prediction.
From figure 7, it is straightforward to see that the effective noise directivity due to the interaction between the instability waves and one shock cell is not of the monopole type, but shows an intrinsic shape that is close to that of the overall screech directivity. This shows that the unique directivity of jet screech is not caused by pure interference between an array of monopoles as assumed by Powell. Therefore, to properly model and understand the directivity of screech, one has to use quantitative models such as the one developed in this paper. However, it is worth noting that this model does not imply that the screech source is localized as a ‘single’ source. What we show is just the effective noise directivity due to interactions between the instability waves and a single shock cell.
To examine the effects of the local growth rate of the instability waves on sound generation, the result obtained using (2.71) with only the real part of 
$\alpha$ considered is shown in figure 8. The SPL is similarly defined by (3.6). Only the fundamental tone and its first harmonic are presented to compare with those reported by Powell (Reference Powell1953a). Considering that the original directivity results reported by Powell (Reference Powell1953a) were presented in the form of schlieren photographs and were therefore not suitable for a direct comparison, only a qualitative comparison is presented. As shown in figure 8, the fundamental tones in all three cases are only reinforced in the upstream direction, while in Powell's experiment sound waves at the fundamental frequency can only be observed propagating upstream, as shown in figure 4 of the original paper (Powell Reference Powell1953a). The model prediction is in agreement with experimental data. In addition, a quick decay also occurs when the observer angle approaches 
$180^\circ$, which is similar to the case when the imaginary part of 
$\alpha$ is not zero. While for the first harmonic, a large lobe perpendicular to the jet flow is predicted by our model. In Powell's experiment, when a reflector was placed, a downstream-propagating sound wave (as shown in figure 5 of the original paper (Powell Reference Powell1953a)) of twice the fundamental frequency emerged. This implied that there was a strong beaming to the side of the jet flow, which is in good qualitative agreement with the prediction.
Figure 8. The directivity of sound in the far field obtained by (2.71). Here 
$r$ is fixed to be 1. Labels (1), (2) and (3) represent the results for a fully expanded jet Mach number of 1.5, 1.3 and 1.2, respectively. Results of (a) the fundamental tone and (b) the first harmonic. The antisymmetric mode of instability waves is taken, and the imaginary part of wavenumber 
$\alpha _{{i}}= 0$. In addition, 
$\mathcal {U}$ in (2.17) is taken to be 1.
To further investigate the effects of the local growth rate of the instability waves on directivity patterns, we change the imaginary part of 
$\alpha$ to 
$2/3$ of its original value. The resulting directivity patterns are shown in figure 9. As can be seen, the fundamental tones in all three cases radiate primarily to the upstream direction, while small lobes appear downstream of the jet flow. Compared with the results shown in figure 7 these lobes are both thinner and weaker, whereas in the case of 
$\alpha _{{i}}=0$ there are no observable lobes downstream of the jet flow, as illustrated in figure 8(a). For the first harmonic, a large lobe appears perpendicular to the jet flow. Compared with the results shown in figure 7, the directivity patterns seem to shrink and move closer to 
$90^\circ$ to the jet. In particular, two small lobes appearing in figure 7(3b) appear to collapse to a single wide lobe, as shown in figure 9(3b).
Figure 9. The directivity of sound in the far field obtained by (2.71). Here 
$r$ is fixed to be 1. Labels (1), (2) and (3) represent the results for a fully expanded jet Mach number of 1.5, 1.3 and 1.2, respectively. Results of (a) the fundamental tone and (b) the first harmonic. The antisymmetric mode of instability waves is taken, and the imaginary part of wavenumber 
$\alpha _{{i}}$ changes to 
$2/3$ of the original value. In addition, 
$\mathcal {U}$ in (2.17) is taken to be 1.
The directivity pattern when a positive imaginary part of 
$\alpha$ is used can be similarly studied. It can be shown that little change occurs when 
$\alpha$ is replaced by 
$\alpha ^{*}$ (
$\alpha ^{*}_{i}=-\alpha _{i}$). So we omit a repetitive discussion for brevity. Figures 7–9 show that the directivity pattern depends on the local growth rate of the instability wave, and a change in the imaginary part of 
$\alpha$ would lead to a corresponding change in the resulting directivity pattern. As mentioned at the beginning of § 3.2, it is not clear exactly where the interaction between instability and shock waves occurs. The amplitude of instability waves may have experienced a growth, a saturation or even decay before reaching the point of interaction. Comparing figures 7–9, we see that the noise directivity is very sensitive to the local growth rate, and this may be used to explain the discrepancies observed across different experiments. For example, Powell's original results showed a clearly dominant radiation only in the upstream direction and a strong 
$90^\circ$ radiation at the first harmonic. This could be explained well if the local instability waves experience a saturation. On the other hand, Walker & Thomas (Reference Walker and Thomas1997), Raman (Reference Raman1999) and Wu, Lele & Jeun (Reference Wu, Lele and Jeun2020) showed that two lobes could be observed at the fundamental frequency, and a relatively weak radiation at 
$90^\circ$ for the first harmonic. This may be explained if the instability waves are in a growth or decay stage at the effective point of interaction.
3.2.2. Directivity patterns from several shock cells
As mentioned at the beginning of § 3.2, some researchers report that screech appears to originate from several shock cells downstream of the jet flow, for example, from the second to the fourth shock cell (Suda et al. Reference Suda, Manning and Kaji1993; Malla & Gutmark Reference Malla and Gutmark2017). Since our model can include multiple shock structures, we can study and compare the sound field produced by the interaction between instability waves and several shock cells with simulations and experiments. Note that to precisely predict the screech amplitude, it is likely that every stage of the feedback loop needs to be considered (Nogueira et al. Reference Nogueira, Jordan, Jaunet, Cavalieri, Towne and Edgington-Mitchell2022b) and the nonlinearity that is inevitable within the loop needs to be included. However, in this paper, we only study the shock-instability interaction and use a linear model (
$\mathcal {U}=1$ in (2.17)). Therefore, the prediction would not be able to match the data in terms of the absolute amplitude. However, we can still plot the predictions and the numerical or experimental data in one figure and focus instead on comparing the shapes of the directivity pattern. The SPL of the model prediction is again defined by (3.6), but rescaled according to the experimental or numerical data. The predictions of the monopole array theory are also included for comparison.
The results are first compared with the study by Wu et al. (Reference Wu, Lele and Jeun2020), where large-eddy simulations were conducted and well validated against the experimental data reported by Alkislar, Krothapalli & Lourenco (Reference Alkislar, Krothapalli and Lourenco2003) and Valentich, Upadhyay & Kumar (Reference Valentich, Upadhyay and Kumar2016). In both the experiment and the numerical simulation, a rectangular nozzle with an aspect ratio of 4:1 was used. The designed Mach number of the nozzle was 1.44, while the Mach number of the fully expanded jet flow was 1.69. It was stated by Wu et al. (Reference Wu, Lele and Jeun2020) that the measured directivity patterns resulted from the interference among spatially distributed sources. Therefore, multiple shock cells are included in our model to facilitate a comparison. Considering that the amplitude of instability waves shows a Gaussian (Freund Reference Freund2011), or more precisely an exponentially modified Gaussian (Maia et al. Reference Maia, Jordan, Cavalieri and Jaunet2019), intensity distribution downstream of the jet flow, in this paper we use three shock cells, the interaction between which and the instability waves leads to different effective source strengths. The middle cell is chosen to have the maximal strength. In front of the middle cell, the instability waves still grow and have not reached the maximum intensity, while after that the instability waves begin to decay but are still of sufficient intensity to generate sound. Thus, the relative strengths of the three interactions are assumed to be 0.45, 1 and 0.7, respectively. Similar assumptions have been also made by Norum (Reference Norum1983) and Berland et al. (Reference Berland, Bogey and Bailly2007). It was known that the effect of varying source strengths on the directivity of the fundamental and the first harmonic was unimportant with regard to the principal lobe, and was appreciable only in the secondary or minor lobes (Kandula Reference Kandula2008). In light of linearity, in what follows we first calculate the sound by one shock cell using our model, and then combine the other two with a spatial phase difference 
$\exp ({\mathrm {i}(\lambda _{0}+\mathrm {real} (\alpha ))2 {\rm \pi}/a})$, where 
$\lambda _{0}=-\omega M_{+}\cos \theta$.
The result is shown in figure 10, where figures 10(a) and 10(b) present directivity patterns of the fundamental tone and the first harmonic, respectively. As can be seen, for the directivity pattern of the fundamental tone, the numerical data show two lobes. The major lobe points to the upstream direction, while the other peaks at around 
$30^\circ$ with a slightly weaker intensity. Note that the radiation intensity seems to decrease as the observer angle approaches 
$0^\circ$. The results from the monopole array theory also show two lobes, which are of equal intensity and peak at 
$180^\circ$ and 
$0^\circ$, respectively. The weaker intensity of the downstream lobe appears not captured in the model. In addition, as 
$\theta$ approaches 
$0^\circ$, the monopole array theory predicts an increasingly large noise radiation, which contradicts the numerical data. On the other hand, we see that the present model predicts a similar major lobe in the upstream direction, and a weaker lobe in the downstream direction. The relative intensity and positions of the two lobes agree better with the numerical data than the monopole array theory. Moreover, the predicted acoustic radiation decreases quickly as the observer angle approaches 
$0^\circ$, which is in good agreement with the numerical data. Note, however, the maximum radiation angle of the downstream lobe appears at around 
$40^\circ$, slightly different from the numerical data. But considering the many assumptions made in the model, such deviation may be deemed acceptable.
Figure 10. Comparison between the numerical results (Wu et al. Reference Wu, Lele and Jeun2020), the present model and the monopole array theory (Powell Reference Powell1953a). The Mach number of the fully expanded jet flow is 1.69. The red solid line denotes the model prediction, the red line with markers the numerical data and the black dashed line the prediction of the monopole array theory. (a) The fundamental tone and (b) its first harmonic.
For the first harmonic, the numerical results exhibit three lobes. Two dominant lobes peak at 
$\theta =30^\circ$ and 
$\theta =83^\circ$, respectively, and one secondary lobe points to 
$155^\circ$. It is evident that a quick decrease occurs as the observer angle approaches 
$180^\circ$ and 
$0^\circ$, respectively. The results obtained by the monopole array theory also show three lobes, which are of the same intensity and peak at 
$180^\circ$, 
$88^\circ$ and 
$0^\circ$, respectively. The corresponding errors compared with the numerical data are 
$+25^\circ$, 
$+5^\circ$ and 
$-30^\circ$, respectively. Moreover, the monopole array theory predicts monotonically increasing acoustic radiations as the observer angle approaches 
$180^\circ$ or 
$0^\circ$, which is not able to match the numerical data. For the present model prediction, a narrow lobe peaks at around 
$84^{\circ }$ with a dominant intensity, and two weaker lobes appear at 
$\theta =31^\circ$ and 
$\theta =135^\circ$, respectively. The corresponding differences compared with the numerical data are 
$+1^\circ$, 
$+1^\circ$ and 
$-20^\circ$, respectively, which is in more satisfactory agreement with the numerical results than those predicted by the monopole array theory. In addition, the rapid decay as the observer angle approaches 
$180^\circ$ and 
$0^\circ$ can be predicted well by this model. Note, however, that the present model cannot correctly predict the relative amplitude of the upstream and downstream lobes. The reason is not yet clear. In summary, although both models are not capable of predicting the screech amplitude, the present model shows a more satisfactory agreement with the numerical data in terms of the peak angle. In addition, unlike the monopole array theory, it can capture the rapid decay of the noise intensity as observer angles approach 
$0^\circ$ and 
$180^\circ$.
The results obtained by the present model are subsequently compared with those of Ponton et al. (Reference Ponton, Manning and Seiner1986), where a series of experiments were conducted at the NASA Langley Research Center. Several microphones were positioned on a circular arc from 
$\theta =30^\circ$ to 
$135^\circ$ at an increment of 
$15^\circ$. The aspect ratio and the designed Mach number of the rectangular nozzle were 3.7 and 1.35, respectively. Two sets of experimental data are chosen for comparison, for which the corresponding fully expanded Mach numbers are 1.5 and 1.6, respectively. Note that the available data only span an observer angle between 
$30^\circ$ and 
$135^\circ$ at an increment of 
$15^\circ$. It is known that very weak noise is radiated in this range, and the peak radiation angles are likely to fall outside this range at the fundamental frequency. Therefore, for a robust comparison we only compare the first harmonic results, where it is known to radiate primarily at side angles. The results predicted by the monopole array theory are also included for comparison. The number and relative strengths of the effective sources are kept the same as those used in figure 10.
As can be seen in figure 11, in the case of 
$M_{-}=1.5$, the experimental data show a major lobe to the side of the jet, and another lobe in the upstream direction with a slightly weaker intensity. In addition, note that in the downstream direction, the acoustic radiation becomes much weaker, as can be seen at an observer angle 
$\theta =30^\circ$. The results from the monopole array theory exhibit three lobes of nearly the same intensity, one in the upstream direction, one in the downstream direction and another to the side of the jet. The agreement with the experimental data is good when 
$60^\circ <\theta <135^\circ$, but much less so in the downstream direction. The present model prediction also exhibits three lobes, one dominant lobe to the side of the jet, one weaker lobe in the upstream direction and another lobe peaking at 
$25^\circ$ with a much weaker intensity. Both the position and the relative intensity of the side and downstream lobes agree well with the experimental data.
Figure 11. Comparison between the experimental data (Ponton, Manning & Seiner Reference Ponton, Manning and Seiner1986), the present model and the monopole array theory (Powell Reference Powell1953a). The fully expanded Mach number is (a) 1.5 and (b) 1.6, while the designed Mach number of the nozzle is 1.35. The red solid line denotes the model prediction, the orange line with markers the experimental results and the black dashed line the prediction of the monopole array theory.
In the case of 
$M_{-}=1.6$, as can be seen in figure 11(b), the experimental data show a dominant lobe to the side of the jet, while another lobe appears slightly downstream with a weaker intensity. Note that a much weaker acoustic wave radiates to the downstream direction. The prediction obtained by the monopole array theory shows three lobes of the same intensity. Again, the agreement with the experimental data is good when 
$60^\circ < \theta < 120^\circ$, but much less satisfactory in the downstream direction. The present model prediction also exhibits three lobes, one dominant lobe to the side of the jet, one weaker lobe in the upstream direction and another peaking at 
$30^\circ$ with a much weaker intensity. Both the position and the relative intensity of the side and downstream lobes agree well with the experimental data. In both the monopole array theory and the present model, it is difficult to draw a conclusion about the agreement when 
$\theta$ approaches 
$135^\circ$ due to the very sparse data points. Because the experimental results only cover an observer angle of 
$30^\circ$ to 
$135^\circ$, the directivity patterns as the observer angle approaches 
$180^\circ$ and 
$0^\circ$ cannot be examined. However, the numerical results in figure 10 show a quick decay as the observer angle approaches 
$180^\circ$ and 
$0^\circ$. A similar phenomenon has been widely reported in numerous experiments for circular nozzles, such as those of Norum (Reference Norum1983) and Powell et al. (Reference Powell, Umeda and Ishii1992). Such an important feature can be well captured by the present model, compared with earlier models.
Note that the present model makes use of a number of linear approximations, for example, both the shock and instability waves are of small magnitude. It is mentioned that supersonic jets may be regarded as weakly imperfectly expanded when 
$|M_{-}^2-M_{1}^2|\le 1$ (Tam & Tanna Reference Tam and Tanna1982; Tam et al. Reference Tam, Jackson and Seiner1985). It can be seen that both the large-eddy simulation and the experiments satisfy this condition. Therefore, the present model may be used to compare with the two cases. However, deviation may occur if intense shocks are involved. In such cases, we would not expect accurate predictions, but it may be possible that some important features of the nonlinear screech may still be captured by the linear model.
In summary, in §§ 3.2.1 and 3.2.2 we show that the effective noise directivity due to the interaction between the instability waves and one shock cell has an intrinsic shape that is close to that of the overall screech directivity, and when incorporating several effective sources, the agreement between the prediction and the experimental data (and numerical results) is more satisfactory. To show the impact of incorporating several sources more clearly, we conduct a quantitative comparison between the directivity patterns obtained using single and multiple shock cells, as can be seen in figure 12. We use the same operating condition as that used in figure 7. Besides, the number and relative strengths of the effective sources are kept the same as those used in figure 11. As shown in figure 12(a), for the fundamental tone, considering multiple shock cells results in two lobes that are both thinner and closer to 
$180^\circ$ and 
$0^\circ$, respectively. This is expected since the three effective sources would lead to constructive interference at 
$180^\circ$ and similar interference may also occur near 
$0^\circ$. For higher harmonics, the upstream and downstream radiations are similarly enhanced. Besides, additional lobes start to appear because of interference. However, we see that the key directivity features from the one-cell interaction remain unchanged. For example, for the first harmonic, the acoustic emission remains strongest to the side of the jet flow, as shown in figure 12(b). Similarly, for the second harmonic, the two dominant lobes are still observable, albeit slightly thinner and tilted. In summary, we see that the final directivity is partly due to the one-cell characteristics and partly due to constructive and destructive interference, as opposed to Powell's assumption that the directivity is solely determined by the interference.
Figure 12. Comparison between the directivity patterns obtained from one effective source and three effective sources, the relative strengths of which are assumed to be (0.45, 1, 0.7). The other operating conditions remain the same as those used in figure 7. Labels (1), (2) and (3) represent the results for a fully expanded jet Mach number of 1.5, 1.3 and 1.2, respectively. Results of (a) the fundamental tone, (b) the first harmonic and (c) the second harmonic. The blue dashed line denotes the directivity pattern from one effective source, while the red solid line represents the results obtained from three sources.
3.3. Near-field pressure and noise generation mechanism
To further examine the noise generation due to the interaction between shock and instability waves, we can calculate the near-field pressure perturbation 
$p_{+}$ by numerically integrating (2.69). In the results shown below, all lengths are non-dimensionalized by the height of the jet 
$D$. The shock spacing can be obtained from (3.1). Figure 13 shows the pressure perturbation immediately outside the jet due to the interaction between instability waves and one shock cell. Figures 13(a)–13(c) present the results at the fundamental frequency and its first and second harmonics, respectively. It can be seen that at the fundamental frequency, the near-field pressure has a dominant distribution in the upstream direction, while at its first harmonic it shows a strong distribution perpendicular to the jet flow. This is in good agreement with the far-field directivity pattern shown in figure 7. For the second harmonic, two major distribution lobes are visible around 
$80^\circ$ and 
$110^\circ$, whereas the radiation at 
$\theta =90^\circ$ is relatively weak. These results are in good agreement with the directivity patterns in the far field, as shown in figure 7(1). When the Mach number of the fully expanded jet changes to 1.3 or 1.2, as shown in figures 14 and 15, respectively, similar agreement between the near field and far field is achieved, and we omit a repetitive description for brevity.
Figure 13. The normalized near-field pressure immediately outside the jet (
$x_{0}$ is the starting position of the effective source) due to one-cell interaction (
$\alpha _{{i}}\neq 0$). The Mach number of the fully expanded jet flow is 1.5. The shock spacing is 2.236; therefore the effective sound source is located between 0 and 2.236. Results at (a) the fundamental frequency, (b) the first harmonic and (c) the second harmonic.
Figure 14. The normalized near-field pressure immediately outside the jet (
$x_{0}$ is the starting position of the effective source) due to one-cell interaction (
$\alpha _{{i}}\neq 0$). The Mach number of the fully expanded jet flow is 1.3. The shock spacing is 1.661; therefore the effective sound source is located between 0 and 1.661. Results at (a) the fundamental frequency, (b) the first harmonic and (c) the second harmonic.
Figure 15. The normalized near-field pressure immediately outside the jet (
$x_{0}$ is the starting position of the effective source) due to one-cell interaction (
$\alpha _{{i}}\neq 0$). The Mach number of the fully expanded jet flow is 1.2. The shock spacing is 1.327; therefore the effective sound source is located between 0 and 1.327. Results at (a) the fundamental frequency, (b) the first harmonic and (c) the second harmonic.
It is interesting to note that for the fundamental frequency the near-field pressure perturbations span the entire shock cell, while for the harmonics they appear to be somewhat localized near the end of the shock, as illustrated in figures 13(b,c), 14(b,c) and 15(b,c). This phenomenon was also reported in a recent experiment conducted by Semlitsch et al. (Reference Semlitsch, Malla, Gutmark and Mihaescu2020). However, since the present model has a periodic nature by construction, we cannot determine whether the sources are physically distributed or localized. But it is interesting to note these near-field behaviours.
Using the model and results of the near-field pressure fluctuations, we now are in a position to examine the noise generation mechanism due to the interaction between shock and instability waves. As illustrated in § 2.4, the velocity potential takes the form
\begin{equation} \phi_{i+}= \frac{U D}{\rm \pi}\int_{-\infty}^{+\infty}D_{1}(\lambda) \exp({-(\mathrm{i}\lambda x+\mathrm{i}\omega t+\gamma_{+}y)})\,\mathrm{d}\lambda\exp({-\mathrm{i}\omega t}). \end{equation}
In the far field, 
$\phi _{i+}$ can be estimated by the saddle point method. The saddle point 
$\lambda _{0}$ is 
$-\omega M_{+}\cos \theta$, where 
$\theta$ represents the observer angle to the downstream direction. It can be shown that 
$D_{1}(\lambda )$ is connected with the Fourier transform of the near-field pressure fluctuations (along a fixed 
$y_{0}$). It is known that in supersonic jets, the phase velocity of the near-field pressure fluctuations along the jet flow can be supersonic relative to the ambient speed of sound, which leads to Mach wave radiation (Tam & Tanna Reference Tam and Tanna1982). As shown in figure 16, the Mach angle satisfies
\begin{equation} \theta^{*}=\arccos\left(\frac{c_{0}}{u_{x}}\right), \end{equation}
where 
$\theta ^{*}$ represents the direction of Mach wave radiation and 
$u_{x}$ and 
$c_{0}$ denote the phase velocities along the jet flow and the radiation direction, respectively. In our case, the phase velocity of the near-field pressure fluctuations in the 
$+x$ direction is equal to 
$-\omega /\lambda$, while the Mach wave has a phase velocity of 
$1/M_{+}$ in the radiation direction. So from (3.8), we can obtain
\begin{equation} \theta^{*}=\arccos \left(-\frac{\lambda}{M_{+}\omega}\right). \end{equation}
It is straightforward to find that 
$\lambda =-\omega M_{+}\cos \theta ^{*}$, which is exactly the same as the saddle point 
$\lambda _{0}$. In fact, the saddle point precisely matches the 
$x$ component of the wavenumber of the Mach wave propagating to the 
$\theta$ direction. As the observer angle (or the Mach wave radiation angle) changes from 
$0$ to 
${\rm \pi}$, the saddle point 
$\lambda _{0}$ changes from 
$-\omega M_{+}$ to 
$\omega M_{+}$ correspondingly. Therefore, the noise radiated at angle 
$\theta$ is directly related to 
$D_{1}(\lambda _{0})$ through the Mach wave mechanism, as shown in (3.7). We may therefore examine the directivities of the sound generation by examining 
$|D_{1}(\lambda )|$ between 
$-\omega M_{+}$ and 
$\omega M_{+}$. Figure 17 shows 
$|D_{1}(\lambda )|$ as a function of 
$\lambda$. We see that sound radiates primarily to the upstream direction (
$\lambda >0$), which is in good agreement with the experimental data.
Figure 16. Schematic of the Mach wave radiation in supersonic jets. The phase velocity of the near-field pressure fluctuations along the jet flow is 
$u_{x}$, while the phase velocity of the radiated sound is 
$c_{0}$.
Figure 17. Plot of 
$|D(\lambda )|$ at the frequency of the fundamental tone. The Mach number of the fully expanded jet flow is 1.3.
However, compared with the result in figure 7, there is no quick decay of 
$D_{1}(\lambda )$ as 
$\lambda \rightarrow \omega M_{+}$, so the sound at 
$180^\circ$ appears to be the strongest. This appears to contradict the results in figure 7. This is because 
$|D_{1}(\lambda _{0})|$ can only show the overall shape of the directivity pattern, but the directivity is in fact determined by 
$|D_{1}(\lambda _{0})\sin \theta |$ instead. We can show this by rewriting (3.7) as
\begin{equation} \phi_{i+}= \frac{U D}{\rm \pi}\int_{-\infty}^{+\infty}D_{1}(\lambda)\exp({-r(\mathrm{i}\lambda \cos\theta+\gamma_{+}\sin\theta)})\,\mathrm{d}\lambda\exp({-\mathrm{i}\omega t}), \end{equation}
where 
$r$ again denotes the radial distance and 
$\theta$ represents the observer angle. We consider the substitution
\begin{equation} \lambda={-}\omega M_{+}\cos(\theta+\beta), \end{equation}with which (3.10) can be rewritten as
\begin{align} \phi_{i+}= \frac{U D}{\rm \pi}\int_{P^{*}}D_{1}(-\omega M_{+}\cos(\theta+\beta))\sin(\theta+\beta)\exp({\mathrm{i}\omega M_{+} r \cos\beta})\,\mathrm{d}\beta\exp({-\mathrm{i}\omega t}), \end{align}
where 
$P^{*}$ denotes the new integration contour in the complex 
$\beta$ plane. We see that the exponent term in (3.12) is independent of 
$\theta$, and now if we use the saddle point method (the saddle point is 
$\beta _{0}=0$), we obtain
\begin{equation} \phi_{i+}= \frac{U D}{\rm \pi}D_{1}(\lambda_{0})\sin\theta F(r,\beta), \end{equation}
where 
$F(r,\beta )$ is independent of 
$\theta$. So the result which determines the far-field directivity pattern is indeed 
$|D_{1}(\lambda _{0})\sin \theta |$.
To better understand this, we take the classical sound field of a monopole as an example. Consider a monopole located at 
$(0,0)$ in the 
$x$–
$y$ plane, as shown in figure 18. Its velocity potential 
$\phi ^{*}$ satisfies the homogeneous Helmholtz equation, i.e.
\begin{equation} \boldsymbol{\nabla}^{2}\phi^{*}+\omega^{2}M_{+}^{2}\phi^{*}=0. \end{equation}
Its axisymmetric solution in the infinite space takes the form 
$\phi ^{*}=H_{0}^{1}(\omega M_{+}r)$, where 
$H^{1}_{0}(\omega M_{+}r)$ is the zeroth-order Hankel function of the first kind and 
$r$ is the radial distance to the sound source. Note that such a solution has a uniform directivity. Although we have obtained the exact solution, we still repeat the calculation process to demonstrate how the factor 
$\sin \theta$ emerges. Similar to (2.69), Fourier transform can be used to solve (3.14), and the solution is
\begin{equation} \phi^{*}(x,y)=\frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty}D^{*}(\lambda)\exp({-(\mathrm{i}\lambda x+\gamma_{+}y)})\,\mathrm{d}\lambda. \end{equation}
Since 
$\phi ^{*}$ is already known, 
$D^{*}(\lambda )$ can be obtained by taking the Fourier transform of 
$H_{0}^{1}(\omega M_{+} r)$ along 
$y=0$. It is straightforward to find that
\begin{equation} D^{*}(\lambda)=\frac{2}{\sqrt{\omega^{2}M_{+}^{2}-\lambda^{2}}}, \end{equation}
where suitable branch of 
$\sqrt {\omega ^{2}M_{+}^{2}-\lambda ^{2}}$ is chosen. In the far field, the saddle point method can be used to estimate (3.15). The saddle point is 
$\lambda _{0}$, and the final result reduces to
\begin{equation} \phi^{*}(r,\theta)=\frac{\sqrt{M_{+}\omega}}{\sqrt{2{\rm \pi}}}D^{*}(\lambda_{0})\sin \theta \frac{\exp({\mathrm{i}\omega(M_{+}r-{\rm \pi}/4)})}{\sqrt{r}}+O(r^{{-}3/2}). \end{equation}
Here 
$D^{*}(\lambda _{0})=2/\omega M_{+}\sin \theta$ and 
$\phi ^{*}(r,\theta )=\sqrt {2/{\rm \pi} \omega M_{+}r}\exp ({\mathrm {i}\omega (M_{+}r-{\rm \pi} /4)})$, which is exactly the far-field approximation of 
$H_{0}^{1}(\omega M_{+} r)$. We can see that the directivity is indeed determined by 
$|D^{*}(\lambda _{0})\sin \theta |$ rather than 
$|D^{*}(\lambda _{0})|$. In fact, the coefficient 
$D^{*}({\lambda }_{0})=2/\omega M_{+}\sin \theta \rightarrow \infty$ as the observer angle 
$\theta \rightarrow 180^\circ$. It is 
$|D^{*}(\lambda _{0})\sin \theta |$ that remains bounded and is independent of 
$\theta$ as expected from the solution 
$H^{1}_{0}(\omega M_{+} r)$.
Figure 18. Schematic of a monopole located at 
$(0,0)$ and its radiated sound.
Figure 17 only shows results for the fundamental frequency at 
$M_{-}=1.3$. Similar results can be obtained for higher harmonics at other Mach numbers. We see that noise is primarily generated through the Mach wave mechanism. Note that Mach wave radiation can also occur in perfectly expanded supersonic jets via jet instability waves. However, in this paper, we focus on the interaction between shock and instability waves, where the near-field fluctuations with supersonic phase speed lead to sound generation. This is not to be confused with the convectional Mach wave radiation due to the jet instability waves in perfectly expanded supersonic jets.
4. Conclusion
An analytical model is developed in this paper to predict the sound arising from the interaction between shock and instability waves in imperfectly expanded two-dimensional jets. Both shock and instability waves are assumed to be of small amplitudes so that linear theories may be used. A vortex-sheet model is used to describe the base jet flow, and two-dimensional Euler equations are subsequently linearized around this base flow to determine the governing equations for shock and instability waves and their interaction, respectively. The interaction between shock and instability waves is determined by solving an inhomogeneous wave equation while simultaneously matching kinematic and dynamic conditions on the vortex sheets. The generated sound in the far field is obtained in a closed form after Fourier transform is used in conjunction with the saddle point method.
The screech frequencies are determined by using the constructive interference assumption proposed by Powell and show good agreement with experimental results. The model can be used to predict the sound due to the instability waves interacting with one shock cell, as well as that with a number of shock cells. The directivity of the sound due to one-cell interaction is shown to resemble that of the total sound field. It is interesting to note that the noise directivity is sensitive to the local growth rate of the instability waves interacting with the shock cells to generate sound and may be used to explain the discrepancies observed across different experiments. When multiple shock cells are included, the present model shows better agreement with experiments and simulations than the monopole array theory. In particular, the present model correctly captures the rapid decay of the acoustic radiation when 
$\theta$ approaches 
$180^\circ$ and 
$0^\circ$. In particular, noise radiation primarily occurs in the upstream direction but becomes weaker as the observer angle gradually approaches 
$180^\circ$, which is in better agreement with experimental results compared with earlier models.
The near-field pressure fluctuation due to the shock–instability interaction is subsequently studied. It is shown that the near-field pressure fluctuation has a distribution that is consistent with the far-field directivity patterns. By examining the wavenumber matching of the near-field pressure, we find that noise is generated primarily through the Mach wave mechanism. It is shown that the model developed in this paper can correctly capture the essential physics and may be used to further study the screech in imperfectly expanded supersonic jets.
Acknowledgements
The second author thanks Professor A.P. Dowling for an earlier stimulating discussion on jet instability waves.
Funding
The authors gratefully acknowledge funding under Marine S&T Fund of Shandong Province for Pilot National Laboratory for Marine Science and Technology (Qiangdao) (no. 2022QNLM010201).
Declaration of interests
The authors report no conflict of interest.
Appendix A
The corresponding velocity, pressure perturbations within the jet and the vortex are
$$\begin{gather} u_{m}=a A \cos(a \beta y)\cos(a x), \end{gather}$$
$$\begin{gather}v_{m}={-}a\beta A \sin(a \beta y)\sin(a x), \end{gather}$$
$$\begin{gather}p_{m}={-}\rho_{-}a U A \cos(a \beta y)\cos(a x). \end{gather}$$
Note that the subscripts in parameters 
$A_{1}$ and 
$a_{1}$ are omitted for clarity. It is worth noting that all these shock-associated perturbations are independent of time. The vortex-sheet deflection at the boundary of jet flow reduces to
\begin{equation} {h}_{m}=\begin{cases} \dfrac{A}{U}\beta\sin\left(\dfrac{1}{2}a\beta\right)\cos(a x), & y=1/2 \\ -\dfrac{A}{U}\beta\sin\left(\dfrac{1}{2}a\beta\right)\cos(a x), & y={-}1/2. \end{cases} \end{equation}Appendix B
The corresponding pressure and velocity are
$$\begin{gather} p_{v}=\mathrm{i}\omega U^{2} \exp({\mathrm{i}(\alpha x-\omega t)})\times \begin{cases} \dfrac{\rho_{+}}{M_{+}^{2}}\,\mathrm{e}^{{-}m_{+}y}, & y>\dfrac{1}{2}\\ \dfrac{\rho_{-}}{M_{-}^{2}}(k_{2}\,\mathrm{e}^{{-}m_{-}y}+k_{1}\,\mathrm{e}^{m_{-}y}), & y\leq \left|\dfrac{1}{2}\right|\\ \dfrac{\rho_{+}}{M_{+}^{2}}k_{3}\,\mathrm{e}^{m_{+}y}, & y<{-}\dfrac{1}{2}, \end{cases} \end{gather}$$
$$\begin{gather}u_{v}=\mathrm{i}\alpha U\exp({\mathrm{i}(\alpha x-\omega t)})\times \begin{cases} \dfrac{1}{M_{+}^{2}}\,\mathrm{e}^{{-}m_{+}y}, & y>\dfrac{1}{2}\\ \dfrac{1}{M_{-}^{2}}\dfrac{\omega}{\omega-\alpha}(k_{2}\,\mathrm{e}^{{-}m_{-}y}+ k_{1}\,\mathrm{e}^{m_{-}y}), & y\leq \left|\dfrac{1}{2}\right|\\ \dfrac{1}{M_{+}^{2}}k_{3}\,\mathrm{e}^{m_{+}y}, & y<{-}\dfrac{1}{2}, \end{cases} \end{gather}$$
$$\begin{gather}v_{v}=U \exp({\mathrm{i}(\alpha x-\omega t)})\times \begin{cases} \dfrac{1}{M_{+}^{2}}m_{+}\,\mathrm{e}^{{-}m_{+}y}, & y>\dfrac{1}{2}\\ \dfrac{1}{M_{-}^{2}}\dfrac{\omega}{\omega-\alpha}m_{-} ({-}k_{2}\,\mathrm{e}^{{-}m_{-}y}+k_{1}\,\mathrm{e}^{m_{-}y}), & y\leq\left|\dfrac{1}{2}\right|\\ \dfrac{1}{M_{+}^{2}}m_{+}k_{3}\,\mathrm{e}^{m_{+}y}, & y<{-}\dfrac{1}{2}. \end{cases} \end{gather}$$The deflection of the jet boundary due to the instability waves is
\begin{equation} {h}_{v}(x,t)=\begin{cases} -\dfrac{\mathrm{i}}{\omega}\dfrac{m_{+}}{M_{+}^{2}}\exp\left({-\dfrac{1}{2}m_{+}}\right) \exp\left({\mathrm{i}(\alpha x- \omega t)}\right), & y=1/2 \\ \dfrac{\mathrm{i}}{\omega}k_{3}\dfrac{m_{+}}{M_{+}^{2}}\exp\left({-\dfrac{1}{2}m_{+}}\right) \exp\left({\mathrm{i}(\alpha x-\omega t)}\right), & y={-}1/2. \end{cases} \end{equation}
