3.1. Building the linear operator

In this study, the global state vector, compactly represented as

(3.1) \begin{align} \boldsymbol{q} = [\xi , u, v, w, p]^T(x, r, \theta , t), \end{align}

includes the specific volume $\xi$ , the velocity components $u$ , $v$ and $w$ in the streamwise, radial and azimuthal directions, respectively, and the pressure $p$ . Specifically, $1/\bar \rho$ , $a_{\infty }$ and $\bar \rho a_{\infty }^2$ are used for the non-dimensionalisation of specific volume, velocities and pressure, respectively, and the perfectly expanded nozzle diameter $D_{\!j}$ is used to non-dimensionalise spatial coordinates. The rationale for the specific-volume/pressure choice introduced in § 2.2 is that the ideal-gas law takes the bilinear form $T = \gamma p \xi$ , which together with the specific-volume substitution $\rho = 1/\xi$ renders every nonlinear term on the right-hand side of the Navier–Stokes equations a product of exactly two components of $\boldsymbol{q}$ ; see Towne (Reference Towne2016) for additional details.

The operator $\boldsymbol{A}$ is constructed by linearising the compressible Navier–Stokes equations around the PIV-derived base flow. Using a normal mode ansatz in the azimuthal and temporal dimensions, the state vector can be expressed as

(3.2) \begin{align} \boldsymbol{q}(x, r, \theta , t) = \hat {\boldsymbol{q}}(x, r)e^{-{i}\omega t + {i}m\theta }, \end{align}

where $\hat {\boldsymbol{q}}(x, r)$ is the spatially varying amplitude of the mode, $m$ is the azimuthal wavenumber and $\omega$ is the angular frequency. Thus, each mode is characterised by a unique pair $(m, \omega )$ , representing its azimuthal and temporal behaviour. In this study, we focus exclusively on the axisymmetric mode ( $m = 0$ ), where the screech phenomenon occurs for the current jet operating conditions.

The computational domain of interest spans $x \times r \in [0, 15] \times [0, 3]$ . The experimental mean flow is extended downstream of the $x = 10$ boundary of the PIV window to $15$ by imposing a self-similar velocity profile from Pope (Reference Pope2000, § 5.2), following prior resolvent and stability analyses of round turbulent jets (Lesshafft et al. Reference Lesshafft, Semeraro, Jaunet, Cavalieri and Jordan2019; Edgington-Mitchell et al. Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021). Sponge layers are applied at the inlet, outlet and far-field boundaries to minimise numerical reflections. The boundary conditions within the sponge region are designed to effectively absorb outgoing waves, thereby reducing reflections and numerical artefacts (Mani Reference Mani2012; Schmidt et al. Reference Schmidt, Towne, Colonius, Cavalieri, Jordan and Brès2017). These sponge layers have been consistently employed in previous numerical studies (Edgington-Mitchell et al. Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021; Prasad & Gaitonde Reference Prasad and Gaitonde2022; Gomez & McKeon Reference Gomez and McKeon2025). The computational domain, including the sponge layers, is discretised with $N_x \times N_r = 700 \times 300$ grid points in the streamwise and radial directions, respectively. A finer grid resolution is adopted in the core region, where smaller flow structures are anticipated, specifically within $x \times r \in [0, 4] \times [0, 0.7]$ , where a non-uniform grid with $350 \times 65$ points is used in the streamwise and radial directions, respectively. Spatial derivatives are computed using fourth-order summation-by-parts finite-difference schemes (Mattsson & Nordström Reference Mattsson and Nordström2004). We verified mesh independence by also computing the leading resolvent modes and gains on a coarser $(N_x \times N_r) = (400, 200)$ grid; differences are small and do not modify the conclusions of the analysis, consistent with the large-scale nature of the dominant KH wavepackets and guided jet modes that are the focus of the present analysis. This grid resolution is also similar to those used in prior global resolvent studies of round jets (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021; Farghadan et al. Reference Farghadan, Martini and Towne2024b ). The linearised Navier–Stokes operator employs a reduced Reynolds number of 1000 to approximate the influence of unmodeled Reynolds stresses (Reynolds & Hussain Reference Reynolds and Hussain1972; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021), consistent with the approach in our earlier work (Farghadan et al. Reference Farghadan, Martini and Towne2024b ). In particular, Pickering et al. (Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021) found $\textit{Re} = 1000$ to be close to the optimal eddy viscosity for a turbulent round jet over a range of frequencies. Moreover, the close agreement between the leading resolvent mode and the leading POD mode at the screech frequency (figure 8) provides an a posteriori check that this choice is reasonable for the present jet.

The norms required for resolvent and harmonic resolvent analyses are defined using Chu’s compressible energy norm (Chu Reference Chu1965) for both forcing and response modes. The weight matrix $\boldsymbol{W} = \boldsymbol{W}_{\!f} = \boldsymbol{W}_{\!q}$ , derived in Appendix A for our specific-volume–pressure formulation, compensates for the non-uniform grid when computing Chu’s energy norm for the non-dimensionalised flow state $\boldsymbol{q}$ . As a result, we define the norm as $|| \boldsymbol{\cdot }|| = || \boldsymbol{\cdot }||_{\!y} = || \boldsymbol{\cdot }||_{\!f}$ , ensuring consistency across the input and output modes. The output and input matrices, $\boldsymbol{C}$ and $\boldsymbol{B}$ , are defined over the domains $x \times r \in [0, 15] \times [0, 3]$ and $x \times r \in [0, 15] \times [0, 1]$ , respectively, excluding the sponge region. These matrices are applied across all analyses, including both resolvent and harmonic resolvent analyses.

3.2. Eigenspectra analysis

We begin our analysis of the screech phenomenon by computing the global eigenspectra of $\boldsymbol{A}$ . Previous studies have shown that the screech mode manifests as a lightly damped, isolated eigenmode in similar jets (Beneddine et al. Reference Beneddine, Mettot and Sipp2015; Edgington-Mitchell et al. Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021).

The eigenvalue spectrum of $\boldsymbol{A}$ for our jet is presented in figure 2 for $\textit{St} \in [0.2, 1]$ . The spectrum consists of several branches and a number of discrete modes. The screech mode, in particular, is expected to manifest at close to the experimentally observed screech frequency. The reported screech frequency is $\textit{St} = 0.714$ , and our analysis identifies the least stable mode (highlighted in red) at $\textit{St} = 0.72$ , corresponding to a relative discrepancy of less than 1 %. In addition to the primary screech mode, three other lightly damped discrete eigenvalues (coloured cyan and green in the figure) are observed; these modes will be analysed shortly.

Figure 2.

Eigenspectrum of the linearised operator around the mean flow. The red eigenvalue is the least-stable mode, which closely matches the observed screech frequency. We show that the two green eigenvalues are additional candidate screech modes created by interactions with different shock-cell wavenumbers. The blue eigenvalue is caused by resonance between upstream- and downstream-travelling acoustic duct modes.

Figure 3 illustrates the real part of the streamwise velocity and pressure components of the eigenmodes corresponding to the four aforementioned discrete eigenvalues. For all subsequent figures of such modes, we show only the real parts, though the modes are inherently complex. For the mode at $\textit{St} = 0.72$ , the pressure field reveals a KH wavepacket localised near the shear layer ( $r \approx 1/2$ ) and a guided jet mode concentrated within the jet core ( $r \lt 1/2$ ), consistent with the interaction between a downstream-travelling KH wavepacket and an upstream-travelling guided jet mode within the screech resonance cycle, as detailed by Edgington-Mitchell et al. (Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021). These characteristics, with KH instabilities dominating in the shear layer and guided modes concentrated within the core, confirm that this mode is a screech mode. The modes at $\textit{St} = 0.66$ and $\textit{St} = 0.77$ exhibit similar features, with wavepackets differing in length and wavenumber due to the change in frequency. We will show in § 3.3 that these discrete modes arise from interactions between KH instabilities and suboptimal wavenumbers of the shock-cell structure (Nogueira et al. Reference Nogueira, Jaunet, Mancinelli, Jordan and Edgington-Mitchell2022a ). Finally, the mode at $\textit{St} = 0.31$ contains both high and low-wavenumber structures in the jet core but lacks a clear KH wavepacket. This structure, as well as the frequency of the eigenvalue, is consistent with resonance between duct-like and guided jet modes, as described by Towne et al. (Reference Towne, Cavalieri, Jordan, Colonius, Schmidt, Jaunet and Brès2017). The localised peak at the centreline near $x \approx 7$ in figure 3(a) reflects the same physics: the contraction of the potential core produces a turning point near that streamwise location, where upstream- and downstream-travelling duct modes resonate (Schmidt et al. Reference Schmidt, Towne, Colonius, Cavalieri, Jordan and Brès2017). From this figure onwards, forcing and response modes are normalised to unit Chu energy and colour bars are accordingly omitted, since the absolute colour values are normalisation-dependent; the sign and spatial structure of each mode (and, for wavenumber spectra of normalised modes, the relative magnitude and peak locations) carry the physical interpretation.

Figure 3.

The eigenmodes associated with isolated eigenvalues of the screeching jet. Panels (a,c,e,g) show the streamwise velocity, and panels (b,d, f,h) show the pressure at frequencies of 0.31 (cyan), 0.66 (green), 0.72 (red) and 0.77 (green), respectively. The colour coding in parentheses corresponds to the spectrum shown in figure 2.

To further investigate the dynamics of these four discrete global modes, we perform a wavenumber transform in the streamwise direction. In what follows, we use $k_{\textit{kh}}$ , $k_{\textit{gj}}$ and $k_s$ to denote the streamwise wavenumbers of the KH wavepacket, the guided jet mode and the shock cells, respectively; their relationship is given formally in (3.3). Figure 4 presents the wavenumber spectra of the pressure component of the eigenmodes at $\textit{St} = 0.31$ , $0.66$ , $0.72$ and $0.77$ (panels a–d), plotting energy as a function of the dimensionless streamwise wavenumbers $k$ and the radial positions $r$ . The white dashed lines represent dimensionless wavenumbers associated with the ambient speed of sound in the upstream and downstream directions. The KH wavepacket and guided jet mode are not periodic in the streamwise direction, instead exhibiting a spatial envelope that leads to energy spreading in wavenumber space. For frequencies near the screeching conditions ( $\textit{St} = 0.66, 0.72, 0.77$ ), the KH wavepacket manifests as a peak observed near the shear layer, with support extending into the core region, at positive wavenumbers in the range $k_{\textit{kh}} \in [5, 8]$ , corresponding to a phase speed of around $0.7 U_{\!j}$ . In contrast, the guided jet mode has its highest amplitude in the core region, appearing at negative wavenumbers $k_{\textit{gj}} \in [-6, -4]$ , which corresponds to waves with a negative phase speed slightly slower than the ambient speed of sound $a_\infty$ (Nogueira et al. Reference Nogueira, Cavalieri, Martini, Towne, Jordan and Edgington-Mitchell2024a ). On the other hand, the mode at $\textit{St} = 0.31$ is confined primarily within the core region of the jet, lacking the shear-layer KH component typical of screech. The wavenumbers are also consistent with those expected for duct modes at this frequency, suggesting that the mode represents a resonance between duct and guided jet modes, as described by Towne et al. (Reference Towne, Cavalieri, Jordan, Colonius, Schmidt, Jaunet and Brès2017). We emphasise that while positive and negative phase speeds do not, in general, imply upstream- and downstream-travelling waves (Briggs Reference Briggs1964; Towne & Colonius Reference Towne and Colonius2015), they do for the wave families involved in screech (Edgington-Mitchell et al. Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021), so we will use this terminology moving forward.

Figure 4.

Wavenumber spectra of the pressure component of the eigenmodes at (a–d) $\textit{St} = 0.31$ , $0.66$ , $0.72$ and $0.77$ , respectively. White dashed lines indicate the acoustic wavenumbers $k_x = \pm \omega$ .

3.3. Resolvent analysis

Next, we compute resolvent modes for the screeching jet. Since the global eigenmode representing the screech was found to be stable, it must be sustained by forcing from background turbulence. Resolvent analysis identifies the optimal forcing and its response, yielding a harmonic screech mode that can be directly compared with data. In particular, the harmonic modes obtained from resolvent analysis provide a better starting point for careful wavenumber analysis than the decaying global eigenmodes (Nichols & Lele Reference Nichols and Lele2011).

Resolvent modes are computed for the frequency range $\textit{St} \in [0.2, 1]$ with a resolution of $ \Delta \textit{St} = 0.01$ using $\boldsymbol{A}$ as the linearised Navier–Stokes operator from § 3.1. Since the focus is on identifying the optimal modes, RSVD-LU is applied with $k = 5$ test vectors and $q=2$ power iterations for all frequencies. The resulting gain spectrum for $m = 0$ is shown in figure 5. The leading gain matches trends observed in prior studies of subsonic jets (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Farghadan et al. Reference Farghadan, Martini and Towne2024b ), while also revealing distinct peaks at $\textit{St} = 0.66$ , $0.72$ and $0.77$ . These prominent peaks are a direct consequence of the interaction between the instability dynamics and the shock-cell structures within the flow, as indicated by the eigenspectrum analysis. The resolvent modes at these frequencies are associated with the screech phenomenon, which constitutes the central focus of this study. The suboptimal gains are more than an order of magnitude lower for most frequencies, and their most prominent feature is a smaller peak at the duct-resonance frequency identified in the global eigenvalue analysis. The suboptimal resolvent response at $\textit{St} = 0.31$ has a structure very similar to the global eigenmode at the same frequency (figure 3), i.e. the duct-resonance mode; the peak appears in the suboptimal gain rather than the leading gain because the duct resonance is essentially uncorrelated with the broadband KH amplification at this frequency and therefore lives in a separate singular direction.

Figure 5.

The resolvent gain $\sigma ^2$ of the screeching jet at $m=0$ , obtained via linearisation about the PIV mean flow: (blue) optimal gain $\sigma _1^2$ ; (orange) suboptimal gains $\sigma _2^2$ and $\sigma _3^2$ .

The streamwise velocity and pressure components of the optimal resolvent responses at these frequencies are shown in figure 6. As expected, these modes exhibit strong similarities to the eigenmodes computed at the same frequencies (Schmid & Henningson Reference Schmid and Henningson2001). A key distinction between the modes at the three frequencies lies in the location of the KH wavepacket envelope, which shifts upstream for higher $\textit{St}$ . This trend is evident in figure 7(a), which plots the absolute value of the pressure along the lip line at $r = 0.5$ . The maximum pressure amplitude occurs progressively closer to the nozzle as $\textit{St}$ increases, with the peak for $\textit{St} = 0.77$ closest to the nozzle, followed by $\textit{St} = 0.72$ , and $\textit{St} = 0.66$ further downstream. This same trend is observed in the eigenmodes at the corresponding frequencies. In all three cases, the presence of a standing wave is observed, as expected given the superposition of upstream and downstream-propagating waves at each frequency. Additionally, figure 7(b–d) displays the wavenumber spectra of the optimal responses, further illustrating the spatial distribution of these modes. The KH instability is clearly observed as a downstream-travelling wave concentrated near the shear layer, while the upstream-travelling guided jet mode is observed closer to the core region. The lowest-frequency mode at $\textit{St} = 0.66$ shows the smallest wavenumber separation between $k_s$ and $k_{\textit{gj}}$ , whereas the higher-frequency cases, $\textit{St} = 0.72$ and $\textit{St} = 0.77$ , exhibit a larger wavenumber separation, suggesting interactions with smaller-scale shock-cell structures located further downstream.

Figure 6.

Resolvent response modes: (a,c,e) show the streamwise velocity, and panels (b,d, f) show the pressure component at frequencies $\textit{St} = 0.66$ , 0.72 and 0.77, respectively.

Figure 7.

Analysis of the resolvent response modes: (a) the absolute pressure along the lip line at $r/D = 0.5$ for $\textit{St} = 0.66$ (red), $\textit{St} = 0.72$ (black) and $\textit{St} = 0.77$ (blue). Panels (b,c,d) show the wavenumber spectra of the optimal pressure responses at frequencies of 0.66, 0.72 and 0.77, respectively. Cyan, magenta and green dashed lines correspond to $k_{\textit{kh}} - k_{s_n}; n =1,2,3$ .

The shock-cell wavenumbers we use below, $k_s = [10.7, 12.3, 13.2]$ for the three leading resolvent modes at $\textit{St} = [0.66, 0.72, 0.77]$ , are obtained from a transform of the standing-wave pattern in the lip-line pressure (figure 7 a), following Edgington-Mitchell et al. (Reference Edgington-Mitchell, Li, Liu, He, Wong, Mackenzie and Nogueira2022); this method circumvents the resolution limit of a direct streamwise transform of the mean pressure (figure 1 b), as we now describe.

As described by Tam & Tanna (Reference Tam and Tanna1982) and Edgington-Mitchell et al. (Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021, Reference Edgington-Mitchell, Li, Liu, He, Wong, Mackenzie and Nogueira2022), the interaction between the streamwise wavenumbers of the KH instability, the shock cells and the guided jet mode is governed by

(3.3) \begin{align} k_{\textit{gj}} = k_{\textit{kh}} - k_s, \end{align}

where $k_{\textit{kh}}$ , $k_s$ and $k_{\textit{gj}}$ represent the wavenumbers of the KH instability, shock cells and guided jet mode, respectively. This equation shows that the downstream-travelling KH wave interacts with the shock cells, and the resulting wavenumber difference gives rise to the upstream-travelling guided jet mode. The rapid variation in shock-cell spacing, as demonstrated in figure 1(a), which shows the mean streamwise velocity profile with weak shock cells and dashed lines indicating their spacing, means that rather than the shock structures being represented by a single value of $k_s$ , they are represented by a broad primary peak and a series of suboptimal peaks. As shown in figure 1(b), the wavenumber spectrum of the mean pressure along the centreline at $r = 0$ exhibits a dominant peak at $k_s \approx 11$ , with cyan dashed lines marking local maxima. The computation of the shock-cell wavenumber $k_s$ , is based on the mean pressure distribution, with a wavenumber resolution of $\Delta k_s = \pi /25 \approx 0.126$ . Due to the gradual reduction in the size of the shock cells along the streamwise direction, the wavelength is not constant but spans a range, leading to a broad peak in the spectrum. This indicates a range of $k_s$ values that can satisfy (3.3). For rapidly varying shock-cell structures, such as those found in jets screeching in the A1 or A2 mode, identifying the first suboptimal peak in the spectrum can be difficult if the primary peak is sufficiently broad. To circumvent the resolution limit of the streamwise transform, we instead identify the appropriate wavenumber through a transform of the standing-wave pattern in the lip-line absolute pressure shown in figure 7; a key result of Edgington-Mitchell et al. (Reference Edgington-Mitchell, Li, Liu, He, Wong, Mackenzie and Nogueira2022) is that the wavenumber of the standing wave always matches the wavenumber of the region of local periodicity that closes the screech loop. The peak wavenumbers of the standing wave for the frequencies $\textit{St} = [0.66,0.72,0.77]$ are $k_s=[10.7,12.3,13.2]$ , whereas the visible peaks of the mean pressure spectra in figure 1(b) are $k_s = [11.0,13.1]$ . Thus, what appears to be the first suboptimal peak in the shock spectra is actually the second suboptimal, with the first suboptimal obscured by the breadth of the primary peak.

Predicted wavenumbers from (3.3) are plotted on figure 7(b,c,d), with the cyan, magenta and green lines corresponding to the wavenumbers $k_{\textit{kh}} - k_{s_n}; n = 1,2,3$ , respectively (the subscript $n$ indexes the primary, first suboptimal and second suboptimal shock-cell wavenumbers in order). It is clear that each of the leading resolvent modes corresponds to a triadic interaction between the KH wavepacket and different regions of quasi-periodicity in the flow; the highest-gain mode at $\textit{St} = 0.72$ , which matches the peak observed in experimental data, results from an interaction between the KH wavepacket and $k_{s_2}$ , while the others correspond to interactions between the primary shock peak and the second suboptimal. To our knowledge, this is the first global stability or resolvent analyses to identify modes associated with multiple shock-cell wavenumbers, and the evidence here rests on the standing-wave extraction described above rather than on a visually resolved three-peak structure in the mean shock spectrum.

In the next section, we will use the resolvent mode at the peak screech frequency as input for a harmonic resolvent analysis. To further confirm that this mode accurately represents screech, we compare the velocity components of the computed resolvent mode and the leading POD mode extracted from experimental data, as shown in figure 8. The close agreement validates the numerical computations and substantiates the identification of the resolvent mode as the screech mode, which serves as the foundation for subsequent analyses. To complete the comparison at the screech frequency, we evaluate the normalised inner product between the eigenmode at $\textit{St} = 0.72$ and the leading resolvent response in the Chu energy norm, $\langle \hat {\boldsymbol{q}}^{{eig}}, \hat {\boldsymbol{q}}_1^R \rangle = 0.998$ , and between the eigenmode and the leading POD mode, $\langle \hat {\boldsymbol{q}}^{ {eig}}, \hat {\boldsymbol{q}}^{\textit{POD}} \rangle = 0.89$ . The first value confirms that the lightly damped eigenmode and the leading resolvent response are collinear at the screech frequency (Schmid & Henningson Reference Schmid and Henningson2001). The second is computed over the streamwise and radial velocity components on the PIV mesh, consistent with the convention used to set $A_s$ , and matches the corresponding resolvent-POD inner product of $0.89$ obtained on the same protocol; the eigenmode and the leading resolvent therefore agree equally well with the POD-derived spatial pattern. Despite this near collinearity, we use the resolvent mode rather than the eigenmode as the periodic base flow in the harmonic resolvent analysis: the eigenmode has a non-zero damping rate and would need to be artificially made periodic for that purpose, whereas the resolvent response at $\textit{St}_s$ is periodic by construction and is therefore the conceptually appropriate starting point.

Figure 8.

Comparison between POD and resolvent modes at the screech frequency: (a,b) POD modes from experimental data and (c,d) optimal resolvent response modes at $\textit{St} = \textit{St}_s$ . Panels (a,c) show the streamwise velocity components, and panels (b,d) display the radial velocity components.

3.4. Harmonic resolvent analysis

We showed in the previous section that resolvent analysis captures interactions between the shock cells, KH waves and guided jet modes. Next, we use the framework of harmonic resolvent analysis to study how the screech mode interacts with other fluctuations in the jet, and, in particular, how it leads to energy transfer to the mean flow and harmonics of the screech frequency, which cannot be captured using resolvent analysis. In what follows, ‘standard resolvent mode’ refers to a mode of the resolvent operator defined in § 2.2, ‘harmonic resolvent mode’ refers to a mode of the harmonic resolvent operator $\boldsymbol{H}$ defined in § 2.3, and ‘optimal’ designates the leading singular direction at a given frequency.

The screech mode oscillates at a well-defined frequency, and a modified mean flow $\bar {\boldsymbol{q}}_{\!T}$ is constructed by superimposing the screech mode $\hat {\boldsymbol{q}}_1^R$ obtained from resolvent analysis onto the original mean flow $\bar {\boldsymbol{q}}$ . This is expressed as

(3.4) \begin{align} \bar {\boldsymbol{q}}_{\!T} = \bar {\boldsymbol{q}} + \boldsymbol{q}_{\!T}, \end{align}

where $\boldsymbol{q}_{\!T} = A_s \mathcal{R}(\hat {\boldsymbol{q}}_1^R e^{{i} 2\pi \textit{St}_s t})$ represents the periodic component. Here, $\textit{St}_s = 0.72$ is the screech frequency, $\mathcal{R}$ denotes the real part of the complex vector and $A_s \approx 0.29$ is the relative amplitude of the screech mode with respect to the mean flow. This relative amplitude was determined by matching the norm of the streamwise and radial velocity components of the resolvent modes to that of the leading POD mode extracted from experimental data shown in figure 8(a,b). Because the POD modes derived from PIV contain only the streamwise and radial velocity components, the matching is restricted to those two components on the PIV mesh; the remaining state components (specific volume, azimuthal velocity, pressure) are not available from PIV and are not included. Using this modified mean flow, a time-periodic linearised Navier–Stokes operator $\boldsymbol{A}_{\!T}$ is constructed. Since the modified mean flow incorporates a single non-zero mode at the screech frequency and we use a bilinear formulation of harmonic resolvent analysis, the resulting $\boldsymbol{T}$ matrix contains only the non-zero components $\hat {A}_{p, 0}$ and $\hat {A}_{p, \pm 1}$ .

The harmonic resolvent modes are computed using the RSVD- $\Delta t$ algorithm (Farghadan et al. Reference Farghadan, Martini and Towne2024b ) with $k=10$ test vectors and $q=1$ power iteration. These parameters enable accurate computation of the first few leading modes. For the time-stepping component, the standard fourth-order Runge-Kutta integration scheme is used with a time step of $\Delta t = 0.011$ , and a transient length of $T_t = 250$ . We also employ the transient-removal strategy proposed by Farghadan et al. (Reference Farghadan, Martini and Towne2024b ), which significantly reduces computational costs. Due to the real-valued nature of $\boldsymbol{A}_{\!T}$ , the negative-frequency modes are the complex conjugates of the positive-frequency modes. As a result, only $\hat {A}_{p, 0}$ and $\hat {A}_{p, 1}$ are retained in memory, and the non-negative fluctuation frequencies are selected as $\varOmega _{\!f} = \{0, \textit{St}_s, 2\textit{St}_s, 3\textit{St}_s, 4\textit{St}_s\}$ . The response is computed up to the $3^{{rd}}$ harmonic $(\textit{St} = 4\textit{St}_s)$ , with convergence observed by the $2^{{nd}}$ harmonic $(\textit{St} = 3\textit{St}_s)$ , indicating that energy contributions from higher harmonics are negligible and do not influence the results. Hence, the analysis is confined to the frequency range within $\varOmega = \varOmega _q = \varOmega _{\!f}$ for the remainder of the paper.

The total CPU time for this analysis is less than 2 hours on 200 cores, using approximately 50 GB of memory. To place the memory requirements in context, the spatial state vector has $N = 5 \times N_x \times N_r = 5 \times 700 \times 300 = 1{\,}050{\,}000$ degrees of freedom per frequency, and the full block-structured harmonic resolvent operator assembled over the nine frequencies $\{-4 \textit{St}_s, -3 \textit{St}_s, \ldots , 3 \textit{St}_s, 4 \textit{St}_s\}$ therefore has $9N \approx 9.45 \times 10^6$ rows. While RSVD-LU achieves comparable CPU efficiency for this problem, its memory requirements are prohibitively high. For the limited frequency range $ \varOmega = \varOmega _q = \varOmega _{\!f} = \{0, \textit{St}_s\}$ , RSVD-LU already consumes around 3.5 TB of memory, making it impractical for broader frequency ranges. In contrast, RSVD- $\Delta t$ offers a more balanced approach, maintaining low memory usage without compromising computational efficiency, making it better suited for the wider frequency range analysed in this study.

In the following paragraphs, we report and discuss the results of the harmonic resolvent analysis. In doing so, it is paramount to keep in mind the meaning of the harmonic resolvent modes: they represent the behaviour of fluctuations (excluding the screech mode) whose dynamics is influenced by interactions with both the mean flow and the screech mode. Because the base flow $\bar {\boldsymbol{q}} + \boldsymbol{q}_{\!T}$ already carries content at $\textit{St} = 0$ and $\textit{St} = \textit{St}_s$ , and the harmonic resolvent fluctuation $\boldsymbol{q}'$ represents motion not already in the base flow, the leading-mode response at these frequencies represents a modification to the assumed base flow. At $\textit{St} = 0$ , this reflects the modification of the mean flow induced by the screech mode interacting with itself. At $\textit{St} = \textit{St}_s$ , the harmonic resolvent response in general need not coincide with the screech mode (other coherent structures can live at the same frequency), and what we find is a structure that resembles the screech mode but with a meaningful phase difference.

Figure 9(a) shows the harmonic resolvent gains. These are plotted as a function of mode index rather than frequency, since each mode may contain all of the frequencies retained in the analysis. The gains exhibit two notable features. First, the leading mode is more than an order of magnitude larger than the second mode; this gain separation suggests that the leading mode is likely to dominate the fluctuation (relative to the mean and the screech mode) for an arbitrary, unknown forcing. Second, despite this gain separation, the leading harmonic resolvent gain is almost an order of magnitude smaller than the leading resolvent gain at the screech frequency (compare with the peak in figure 5). This is consistent with the fact that the screech mode is by far the largest fluctuation (relative to the mean) observed in the experimental data.

Figure 9.

Harmonic resolvent gains and energy distribution for the screeching jet: (a) the six leading harmonic resolvent gains. For comparison, the horizontal dashed line shows the gain of the leading resolvent mode. (b) The energy spectrum of the optimal harmonic resolvent response (red) and forcing (blue) modes. The $q$ -norm of the response (with all frequencies combined) is 1.

The energy distribution of the optimal harmonic resolvent mode as a function of frequency, shown in figure 9(b), highlights the role of the screech mode in distributing energy across frequencies. Unlike the optimal resolvent mode, which oscillates strictly at the screech frequency, the harmonic resolvent mode contains significant energy content at the screech frequency, while also feeding energy into the zeroth and higher harmonics through triadic interactions with the screech mode. In particular, the zeroth and first harmonics contain nearly half of the total energy of the mode. This observation underscores the pivotal role of the screech mode in redistributing energy across the frequency spectrum.

Next, we analyse the spatial structures of the leading harmonic resolvent mode at each frequency and compare them with those of the standard resolvent modes at the same frequencies to understand their physical characteristics. Recall that, due to the large gain separation in the harmonic resolvent gain spectrum, the leading harmonic resolvent mode shapes represent the dominant physical structures that we expect to observe in the fluctuation field (relative to the mean flow and screech mode). The standard resolvent modes were computed using the RVSD-LU algorithm with $k=5$ test vectors and $q=2$ power iterations. Results are shown in figure 10 and are discussed in order of decreasing energy content within the harmonic resolvent mode.

Figure 10.

Comparison between optimal resolvent and harmonic resolvent response modes: (a,c,e,g) resolvent pressure modes and (b,d, f,h) harmonic resolvent pressure modes at frequencies of $\textit{St} = 0, \textit{St}_s, 2\textit{St}_s$ and $3\textit{St}_s$ , respectively.

At the screech frequency, $\textit{St} = \textit{St}_s$ , the resolvent mode features the signatures of a downstream-travelling KH wavepacket and an upstream-travelling guided jet mode, as discussed in detail in § 3.3. The harmonic resolvent mode is visually similar, exhibiting the same two waves. A quantitative comparison confirms this close agreement; projecting the two modes onto one another yields a normalised inner product of 0.97, indicating that they are nearly collinear in the chosen norm. Together, these results indicate that the easiest structure to force within the screeching jet is a modification to the screech mode itself. Physically, this fluctuation is enabled by the same resonance loop underlying the screech mode via its interaction with the mean flow, but is at the same time modulated by the presence of the screech mode, leading to energy redistribution to other frequencies and a lower gain than the screech mode itself.

At zero frequency, $\textit{St}=0$ , the resolvent mode contains modifications to the shock cells throughout the potential core as well as likely unphysical high-wavenumber noise within the downstream turbulent region of the jet. The gain of this mode is low, making it both difficult to converge numerically and unlikely to play a prominent role in flow physics. In contrast, the harmonic resolvent mode consists entirely of a modification to the shock cells focused in the region where the screech mode is active, and the amplitude of this component of the mode is significant.

At the first harmonic, $\textit{St} = 2 \textit{St}_s$ , the resolvent mode consists of a KH wavepacket near the nozzle (consistent with the higher frequency) and guided jet modes with higher radial order (Towne et al. Reference Towne, Cavalieri, Jordan, Colonius, Schmidt, Jaunet and Brès2017). The harmonic resolvent mode is completely different, consisting of a wavepacket in the region where the screech mode is active and acoustic waves propagating perpendicular to the jet axis. These acoustic waves have been observed in experiments (Powell Reference Powell1953b ) but have not been previously captured within linearised Navier–Stokes analyses. The present results show that the interaction of fluctuations with the screech mode can generate the observed acoustic signature. At the second harmonic, $\textit{St} = 3 \textit{St}_s$ , the resolvent mode consists of direct acoustic radiation. Again, the harmonic resolvent mode is completely different, consisting of a concentrated wavepacket in the region where the screech mode is active and a beam of acoustic radiation that is again consistent with experimental observations (Tam et al. Reference Tam, Parrish and Viswanathan2014). Overall, the differences between resolvent and harmonic resolvent modes indicate the consequential impact of the screech mode on the dynamics of other fluctuations in the jet, including the redistribution of energy to other frequencies.

Figure 11 compares the optimal forcing modes from resolvent and harmonic resolvent analyses. At the screech frequency, $\textit{St} = \textit{St}_s$ , the forcing modes are nearly identical, consisting of slanted structures in the near-nozzle region that effectively excite the KH waves (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018) and more distributed acoustic-like structures that effectively excite the guided jet mode due to the self-adjoint nature of acoustic waves. The similarity between the resolvent and harmonic resolvent forcing mode reiterates the point that the best way to generate fluctuations at the screech frequency is to excite the screech resonance, even with the redistribution of energy due to the presence of the screech mode itself.

Figure 11.

Comparison between optimal resolvent and harmonic resolvent forcing modes: (a,c,e,g) resolvent pressure-forcing modes and (b,d, f,h) harmonic pressure-forcing modes at frequencies of $\textit{St} = 0, \textit{St}_s, 2\textit{St}_s$ and $3\textit{St}_s$ , respectively.

The forcing modes at zero frequency, $\textit{St}=0$ , are consistent with the response modes at the same frequency, with the former spread throughout the potential core and the latter concentrated in the region where the screech mode is active. At both of the harmonic frequencies, $\textit{St} = 2 \textit{St}_s$ and $3 \textit{St}_s$ , the resolvent forcing modes are dominated by the direct forcing of acoustic waves. In contrast, the tilted structures that effectively force the KH waves are more prominent in the harmonic resolvent forcing modes. Given the low gains of the resolvent modes at the non-screech frequencies and the differences in the forcing modes, we conclude that the harmonic resolvent response modes at non-screech frequencies are primarily driven by energy transfer due to the presence of the screech mode, rather than direct linear amplification.

Finally, the wavenumber spectra of the optimal harmonic resolvent response modes, shown in figure 12, provide insights into the wave interactions responsible for these modes. At the screech frequency, $\textit{St} = \textit{St}_s$ , the wavenumber spectrum is qualitatively similar to that of the resolvent screech mode (shown previously in figure 7); there are minor differences in the peak wavenumber of the KH wavepacket, but little else distinguishes the two modes. Weaker peaks appear at the negative wavenumber $k_{\textit{kh}}-(k_{s_{1}}+k_{s_{2}})$ and the positive wavenumber $k_{\textit{kh}} + k_{s_{1}}$ that are confined within the core of the jet, consistent with downstream-travelling duct modes (Towne et al. Reference Towne, Cavalieri, Jordan, Colonius, Schmidt, Jaunet and Brès2017).

Figure 12.

The wavenumber spectra of the pressure component of the optimal harmonic resolvent response at (a–d) the frequencies $\textit{St} = 0, \textit{St}_s, 2\textit{St}_s$ and $3\textit{St}_s$ , respectively. The vertical dashed lines indicate wavenumbers for (white) acoustic waves and (cyan, magenta, green) $k_{\textit{kh}} - k_{s_n}$ ; $n=1, 2, 3$ .

At $\textit{St}=0$ , the harmonic resolvent mode is real valued and mirror symmetric about zero, representing the stationary modification of the mean flow and shock cells by the unsteady structures active at the screech frequency. The wavenumber spectrum displays dominant peaks at shock-associated wavenumbers (e.g. $k_{s_{2}}$ and $k_{s_{3}}$ ) as well as high-amplitude peaks at various sum and difference combinations. While the individual peaks suggest complex interactions, the spectrum exhibits a consistent spacing of approximately $\Delta k = 11$ , beginning at $k_{\textit{mod}} = 15$ . This starting wavenumber corresponds to the self-interaction of the KH wavepacket observed at the screech frequency ( $2k_{\textit{kh}}$ , where $k_{\textit{kh}} \approx 7.5$ ). This spectral signature confirms a two-way coupling mechanism: the shock structures modulate the KH wavepacket (Nogueira et al. Reference Nogueira, Self, Towne and Edgington-Mitchell2022c ), and the nonlinear self-interaction of this modulated wavepacket in turn modifies the mean flow.

At the first harmonic, $\textit{St} = 2\textit{St}_s$ , in addition to the KH wavepacket appearing at twice the wavenumber of the fundamental, there are peaks associated with both $k_{\textit{kh}}-k_{s_{1}}$ and $k_{\textit{kh}}-k_{s_{2}}$ ; the former has strong support outside of the jet, whereas the latter has the highest amplitude in the jet core. These wavenumbers, close to zero, are consistent with the directivity of the beams of acoustic radiation observed in figure 10(f) and in Tam et al. (Reference Tam, Parrish and Viswanathan2014). Interestingly, the results here suggest that, while $k_{s_{2}}$ is responsible for closing the resonance loop, much of the sideline radiation of the harmonic is in fact associated with an interaction between $k_{\textit{kh}}$ and $k_{s_{1}}$ . There is also a component of upstream radiation associated with $k_{\textit{kh}}-(k_{s_{1}}+k_{s_{2}})$ , which we again suggest is the signature of a shock-modulated wavepacket undergoing triadic interaction with another shock wavenumber.

At the second harmonic, $\textit{St} = 3\textit{St}_s$ , the KH wavepacket peaks at $k_x \approx 21$ , approximately triple the value observed at $\textit{St} = \textit{St}_s$ , as expected for a second harmonic. Here, the peaks associated with the interactions $k_{\textit{kh}}-k_{s_{1}}$ and $k_{\textit{kh}}-k_{s_{2}}$ are both relatively weak; there is some redistribution of energy to these wavenumbers, but no clear local peak. The stronger peaks in fact occur at $k_{\textit{kh}}-(k_{s_{1}}+k_{s_{2}})$ and $k_{\textit{kh}}-(2k_{s_{1}}+k_{s_{2}})$ . We suggest that these results can be interpreted as the KH wavepacket, modulated by its interaction with $k_{s_{2}}$ , interacting with $k_{s_{1}}$ and its harmonic, respectively.

3.5. Nonlinear forcing

The harmonic resolvent analysis in the previous section examined the response of the jet to the optimal forcing, which is the sensible approach when the true nonlinear forcing $\hat {\kern-3pt \boldsymbol{f}}$ is unknown. However, we showed in § 2.4 that for a bilinear system, the forcing is not entirely unknown. Specifically, our bilinear formulation allows us to directly study how the nonlinear self-interaction of the screech mode, $\hat {\kern-3pt \boldsymbol{f}}_{nl} = \boldsymbol{\mathcal{B}}(\boldsymbol{q}_{\!T}, \boldsymbol{q}_{\!T})$ , impacts the other fluctuations in the jet. This approach helps elucidate the role of the screech mode in redistributing energy to harmonic frequencies. Throughout this section, the action of $\boldsymbol{H}$ on $\hat {\kern-3pt \boldsymbol{f}}_{nl}$ is computed using time-stepping analogous to the RSVD- $\Delta t$ algorithm.

According to (2.28), the forcing term $\hat {\kern-3pt \boldsymbol{f}}_{nl}$ is non-zero at $\textit{St} = 0$ and $\textit{St} = \pm 2\textit{St}_s$ due to the concentration of $\hat {\boldsymbol{q}}_1^R$ at $\textit{St} = \pm \textit{St}_s$ and the non-zero nature of $\hat {\boldsymbol{A}}_i$ for $|i| \leqslant 1$ . Figure 13(a) illustrates this behaviour by showing the energy spectrum of both $\hat {\boldsymbol{q}}_1^R$ and $\hat {\kern-3pt \boldsymbol{f}}_{nl}$ across different frequencies. The normalisation $||\hat {\boldsymbol{q}}_{\pm 1}||^2 = 1$ ensures that the total energy of $\hat {\boldsymbol{q}}_1^R$ is $||\hat {\boldsymbol{q}}_1^R||^2 = 2$ . More generally, the total energy of any state $\boldsymbol{q}$ across all frequencies is defined as

(3.5) \begin{align} || \hat {\boldsymbol{q}} || ^ 2 = \sum _{i = -3}^{3} || \hat {\boldsymbol{q}}_i ||^2, \end{align}

where $\hat {\boldsymbol{q}}_i$ represents the state at the ith harmonic of the screech frequency. The total energy of the nonlinear forcing term is computed $||\hat {\kern-3pt \boldsymbol{f}}_{nl}||^2 = 38.0$ . Figure 13(b) presents the energy distribution of the nonlinear response $\hat {\boldsymbol{q}}_{nl}$ , revealing that its total energy is $||\hat {\boldsymbol{q}}_{nl}||^2 = 25.6$ , which is lower than the input forcing norm. This indicates that the forcing undergoes some degree of attenuation. Notably, while the nonlinear forcing is zero at $\textit{St} = \textit{St}_1$ and $\textit{St}_3$ , the triadic interactions retained within harmonic resolvent analysis redistribute energy to these frequencies.

Figure 13.

Energy distribution of quantities in the nonlinear analysis: (a) $\hat {\boldsymbol{q}}_1^R$ (black), $\hat {\kern-3pt \boldsymbol{f}}_{nl}$ (blue), (b) $\hat {\boldsymbol{q}}_{nl}$ (red).

The pressure modes of both $\hat {\kern-3pt \boldsymbol{f}}_{nl}$ and $\hat {\boldsymbol{q}}_{nl}$ are depicted in figure 14. While the nonlinear forcing modes at $\textit{St} = 0$ and $\textit{St} = \pm 2\textit{St}_s$ differ significantly from the optimal harmonic forcing modes in figure 11, the pressure response modes of $\hat {\boldsymbol{q}}_1^H$ (figure 10) and $\hat {\boldsymbol{q}}_{nl}$ closely resemble each other. To quantify this agreement, we compute the inner products between $\hat {\boldsymbol{q}}_{1, i}^H$ and $\hat {\boldsymbol{q}}_{nl, i}$ for $|i| \leqslant 3$ , obtaining values greater than 0.95 in all cases, indicating a strong alignment at each individual frequency. However, the inner product $\langle \hat {\boldsymbol{q}}_1^H, \hat {\boldsymbol{q}}_{nl} \rangle = 0.36$ suggests a phase mismatch at some frequencies. Here, the values greater than $0.95$ are normalised inner products $\langle \hat {\boldsymbol{q}}_{1,i}^H, \hat {\boldsymbol{q}}_{nl,i} \rangle$ evaluated frequency by frequency, while the value $0.36$ is the normalised inner product of the full frequency-concatenated modes. The two differ because the harmonic resolvent SVD selects the relative phase between frequency blocks that maximises the total response, whereas $\hat {\boldsymbol{q}}_{nl}$ inherits whatever phase relationship the known nonlinear self-interaction produces; this phase difference is meaningful, not arbitrary. The two analyses are therefore complementary: the harmonic resolvent SVD identifies the optimal receptive direction $\hat {\boldsymbol{q}}_1^H$ independently of the forcing, while the present analysis reveals the response to the specific known nonlinear self-interaction, which is best interpreted as a correction to that optimal direction rather than as exciting it directly.

Figure 14.

Spatial structure of quantities in the nonlinear analysis: (a,c,e,g) forcing pressure modes and (b,d, f,h) corresponding response pressure modes at frequencies $\textit{St} = 0, \textit{St}_s, 2\textit{St}_s$ and $3\textit{St}_s$ , respectively.

To summarise, this analysis shows that the forcing provided by the direct nonlinear self-interaction of the screech mode, along with its triadic interactions with other frequencies, is sufficient to explain the redistribution of energy to other frequencies and the associated properties of these fluctuations. These include modifications to the shock cells and screech mode, and the experimentally observed acoustic beams at the harmonic frequencies. No appeal to an unknown forcing by background turbulence is necessary to explain or model these phenomena.