Resolvent analysis has demonstrated encouraging results for modeling coherent structures in jets when compared against their data-educed counterparts from high-fidelity large-eddy simulations (LES). We formulate resolvent analysis as an acoustic analogy that relates the near-field resolvent forcing to the near- and far-field pressure. We use an LES database of round, isothermal, Mach 0.9 and 1.5 jets to produce an ensemble of realizations for the acoustic field that we project onto a limited set of resolvent modes. In the near-field, we perform projections on a restricted acoustic output domain, $r/D=[5,6]$, while the far-field projections are performed on a Kirchhoff surface comprising a 100-diameter arc centered at the nozzle. This allows the LES realizations to be expressed in the resolvent basis *via* a data-deduced, low-rank, cross-spectral density matrix. We find that a single resolvent mode reconstructs the most energetic regions of the acoustic field across Strouhal numbers, $St=[0\u22121]$, and azimuthal wavenumbers, $m=[0,2]$. Finally, we present a simple function that results in a rank-1 resolvent model agreeing within 2 dB of the peak noise for both jets.

## I. INTRODUCTION

The goal of this work is to develop jet-noise models founded upon the physics of turbulent flows that are both low-rank and that provide insights into the mechanisms primarily responsible for noise generation. Resolvent analysis (McKeon and Sharma, 2010), also known as input-output analysis (Jovanović, 2021), provides a useful framework for achieving these goals. The central idea of the resolvent framework is similar to that of an acoustic analogy (Goldstein, 2003; Lighthill, 1952), whereby a forcing term, related to the statistics of the hydrodynamic near-field turbulence, gives rise, through a linear operator, to the observed far-field sound. The resolvent framework differs in two important ways. First, the operator is decomposed into its singular components that represent the maximal amplification between the forcing and the output. This permits the resulting acoustic field to be described as low rank, and thus limits the number of forcing statistics that must be modeled. Second, the full linearized Navier-Stokes equations are used as the propagator, and we seek a modal basis that represents *both* near and far-field coherent structures.

Before recent advances in computational power, the idea of modeling both the hydrodynamic component along with the acoustics would have been seen as both unnecessary and computationally taxing. However, the ability to resolve both components of the flow is in fact a benefit. Starting with the experimental findings of Mollo-Christensen (1967) and Crow and Champagne (1971), it has become clear that coherent structures in the hydrodynamic near-field are directly responsible for far-field sound (Jordan and Colonius, 2013). These structures take the spatiotemporal form of wavepackets and have been found to be the dominant source for aft-angle sound (Jordan and Colonius, 2013), as well as partial contributors to sideline noise (Jeun and Nichols, 2018; Papamoschou, 2018). These wavepackets may be linked to the early works of Crighton and Gaster (1976) (and Michalke, 1977), who hypothesized that coherent structures could be described as linear instability modes of the mean flow *via* modal analysis. However, it has now become apparent that the correct representation of wavepackets is that of a highly-amplified response to turbulent fluctuations, which is directly found *via* the resolvent framework.

Resolvent analysis uses the Singular Value Decomposition (SVD) to decompose the linear resolvent operator, identifying sets of orthogonal *forcing/input* and *response/output* modes, and ranking them in terms of the corresponding energetic gain between the forcing and response. This is particularly important as it allows our model to self-select the most relevant amplification mechanisms for noise generation. This allows for a natural truncation of the resolvent basis that produces a reduced-order model, or in other words, a reduced-rank acoustic analogy.

Several studies have applied resolvent analysis to develop low-rank jet models (Cavalieri *et al.*, 2019; Jeun *et al.*, 2016; Lesshafft *et al.*, 2019). The existence of relatively low-rank responses in round, turbulent jets was shown by Schmidt *et al.* (2018), with significant agreement between structures found through spectral proper orthogonal decomposition (SPOD) (Towne *et al.*, 2018) of a high-fidelity experimentally-verified large-eddy simulations (LES) of jets (Brès *et al.*, 2017; Brès *et al.*, 2018). Of particular relevance to this study are “acoustic resolvent modes” induced by performing resolvent analysis with an output domain defined over a region where fluctuations are purely acoustic. Through implementation of an acoustic output domain, resolvent analysis is able to filter out energetic, but acoustically irrelevant structures in the near-field. Jeun *et al.* (2016) performed such an analysis and found that for a Mach 1.5 jet, at Strouhal number *St =* 0.33 and azimuthal wavenumber *m* = 0, the first resolvent mode reconstructs 57% of the acoustic energy, but through inclusion of the next 23 resolvent modes, the reconstruction improved to 70% of the acoustic energy. This study looks to perform a similar analysis, in that we compute many acoustic resolvent modes and assess how well they reconstruct the acoustic energy. However, we also look to reduce the rank of the far-field significantly with the use of an eddy-viscosity model (Pickering *et al.*, 2021) and generalize the performance of the resolvent framework across frequencies $St=0\u22121$, azimuthal wavenumbers $m=[0\u22122]$, and for two turbulent jets at Mach numbers of 0.9 and 1.5.

For a resolvent jet model to fully reconstruct flow statistics, and in this case those of the acoustic field, a resolvent-based model must incorporate sub-optimal modes (Schmidt *et al.*, 2018) and correctly describe correlations (i.e., covariance) between modes inherent to turbulent flow (Towne *et al.*, 2020). These correlations are analogous to the concept of “jittering” used to describe temporal modulations of acoustic sources, which has been shown to be critical for accurately describing the acoustic field in turbulent jets (Cavalieri *et al.*, 2011). In our approach, such temporal modulations, or jittering, may be represented through second-order statistics *via* the statistical representation of the resolvent operator (Towne *et al.*, 2018)

where $Syy$ and $Sff$ are the cross-spectral density (CSD) tensors of the response and the forcing, respectively, and ** R** is the resolvent operator. This equation shows that if the forcing CSD, describing spatial correlations, can be modeled (Towne

*et al.*, 2017; Zare

*et al.*, 2017), then the resolvent operator identically reconstructs the flow statistics, $Syy$. If the forcing were spatially uncorrelated, $Sff=\Lambda $, where $\Lambda $ is a diagonal matrix, then the eigenvectors of $Syy$, which are the SPOD modes of the outputs, are aligned with the eigenvectors of $RR*$ (Towne

*et al.*, 2018), or the response modes of the resolvent operator,

**. However, the uncorrelated condition is rarely met, resulting in discrepancies between resolvent and SPOD modes that must be resolved through modeling $Sff$.**

*R*One approach for modeling $Sff$ has been through the inclusion of a turbulence model within the resolvent operator. This approach has been implemented *via* an eddy-viscosity model in several flow configurations, from wall-bounded (Hwang and Cossu, 2010; Morra *et al.*, 2019) to free shear flows (Pickering *et al.*, 2021). The latter study, quantifying the effect on turbulent jet modeling, found that the use of an eddy-viscosity model [utilizing only quantities available from Reynolds-averaged Navier-Stokes (RANS) models] significantly improved the agreement between SPOD and resolvent modes, thus reducing the effort required to model the effective $Sff$ by diminishing the magnitude of the off diagonal terms. We utilize the same eddy-viscosity model in the present work to better model the acoustic field.

This paper explores an approach to describe the coupling between resolvent modes that is necessary for reconstructing the acoustic field with a minimal set of resolvent modes. The coupling provides directional and energetic variability in acoustic radiation inherently important for noise prediction (Cavalieri *et al.*, 2011). Determination of the coupling between modes is performed by leveraging an ensemble of LES realizations which are projected on to a limited (i.e., low-rank) set of acoustic resolvent modes. From these projections, we attain a (drastically) reduced-order cross-spectral density between the retained modes—a Hermitian, frequency-dependent matrix of size *n* ×* n*, where *n* denotes the number of retained modes, that accurately represents the acoustic field.

Organization of the manuscript is as follows. We first briefly describe the LES databases used, the main details pertaining to resolvent analysis, and present the statistical description of the resolvent framework for reconstructing the acoustic field and estimating the reduced order covariance matrix in Sec. II. In Sec. III, we present resolvent modes and LES reconstructions in the resolvent basis for one frequency-wavenumber pair for the Mach 1.5 jet before generalizing the approach to both jets over $St=[0,1]$ and $m=[0,2]$, and to both the near- and far-field acoustic regions. In the near-field section we compare the impact of including a RANS eddy-viscosity model to the resolvent operator and find it presents a significantly more efficient resolvent basis. We then present results for the far-field, along an arc at 100*D* from the nozzle, and show that reconstructions for both jets may be found using only the optimal resolvent mode. Finally, we conclude with a discussion on how the correct forcing coefficients may be estimated for a predictive jet noise model.

## II. METHODS

### A. Large Eddy simulation database

The LES database and resolvent analysis are fully described in Schmidt *et al.* (2018) and Towne *et al.* (2018). Transonic (Mach 0.9) and supersonic (Mach 1.5) jets were computed using the flow solver “Charles”; details on numerical methods, meshing, and subgrid-models can be found in Brès *et al.* (2018) and Brès *et al.* (2017) along with validation cases conducted at PPRIME Institute, Poitiers, France for the Mach 0.9 jet (Brès *et al.*, 2018). The Mach 0.9 and 1.5 jets have Reynolds numbers of $Rej=\rho jUjD/\mu j=1.01\xd7106$ and $Rej=1.76\xd7106$, respectively, where subscript *j* gives the value at the center of the jet, *ρ* is density, *μ* is viscosity, and *M _{j}* is the Mach number $Mj=Uj/cj$, with

*c*as the speed of sound at the nozzle centerline.

_{j}Throughout the manuscript, variables are non-dimensionalized by the mean jet velocity *U _{j}*, jet diameter

*D*, and pressure $\rho jUj2$, with the resulting equation of state $p=\rho T/\gamma Mj2$, with

*T*denoting temperature and

*γ*the ratio of specific heats. Frequencies are reported in Strouhal number, $St=fD/Uj$, where

*f*is the frequency in Hertz. The database consists of 10 000 snapshots separated by $\Delta tc\u221e/D=0.2$ and 0.1 for the $Mj=0.9$ and $Mj=1.5$ jets, respectively, with $c\u221e$ as the ambient speed of sound, and interpolated onto a structured cylindrical grid $x,r,\theta \u2208[0,30]\xd7[0,6]\xd7[0,2\pi ]$, where

*x*,

*r*,

*θ*are streamwise, radial, and azimuthal coordinates, respectively. Variables are reported by the vector

where *u _{x}*,

*u*, $u\theta $ are the three cylindrical velocity components.

_{r}To generate an ensemble of flow realizations for computing statistical averages, the LES database of 10 000 snapshots is segmented into bins of 256 snapshots, with an overlap of 75%, and under the implementation of a Hamming window, resulting in 153 realizations of the flow. Each realization is then decomposed in the azimuthal direction and in time. The temporal decomposition provides a resolution of *St =* 0.026 and *St =* 0.0217 the $Mj=1.5$ and $Mj=0.9$ jets, respectively, and the azimuthal decomposition is valid up to *m =* 68; however, the acoustically relevant azimuthal wavenumbers are much smaller (Juve *et al.*, 1979) and only azimuthal wavenumbers $m=[0\u22122]$ are considered in this paper.

Considering the LES database only extends to $r/D=6$, we implement a Kirchhoff surface (as described in Freund, 2001), to the azimuthally and temporally transformed realizations of the flow. In doing so, we create an ensemble of far-field realizations located along an arc, with angle $\varphi $, of 100*D* from the nozzle at each frequency and azimuthal wavenumber. As done in Brès *et al.* (2017), and associated experiments (Schlinker *et al.*, 2009; Schlinker *et al.*, 2008), we specifically compute the acoustics for the aft-angle sound from $\varphi =100\u2212160$ and find our acoustic far-field is in close agreement (within 2 dB) with the far-field of the LES calculation.

### B. Resolvent analysis

For the round, statistically-stationary, turbulent jets considered in this manuscript, the compressible Navier-Stokes, energy, and continuity equations are linearized *via* a standard Reynolds decomposition (see Pickering *et al.*, 2021, Appendix B for a detailed description of the governing equations in cylindrical coordinates) and Fourier transformed both in time and azimuthally to the compact expression

where $\omega =2\pi St$ is the frequency, *m* is the azimuthal wavenumber, **I** is the identity matrix, $Am$ is the frequency independent linear operator, $Lm,\omega $ is the total forward linear operator, $qm,\omega $ is the response in each variable, and $nm,\omega $ constitutes the nonlinear forcing. Mean-flow quantities used in the operator are derived from a RANS model, fitted closely to the LES mean flow. Although the mean flows are similar, the computation of a RANS model, using the standard $\kappa \u2212\u03f5$ closure equations, also provides an eddy-viscosity field that may be included in the resolvent operator. This is done following the results of Pickering *et al.* (2021) that presented substantially improved agreement between SPOD and resolvent modes with the inclusion of an eddy-viscosity model. The eddy-viscosity used here is computed as $\mu T=cC\mu k2/\u03f5$, where *c* and $C\mu $ are scaling constants (*c =* 0.2, $C\mu =0.0623$ for the $Mj=0.9$ and $C\mu =0.0554$ for $Mj=1.5$ jet), *k* is the turbulent kinetic energy field, and *ϵ* is the turbulent dissipation field. We stress that the values of $C\mu $ used here were selected to reproduce the LES mean flow and to allow for a demonstration of the ability of quantities available in RANS to be used to compute accurate resolvent modes. These values should not be considered general or as recommended values of $C\mu $ for general purposes.

Continuing with the derivation of the resolvent/input-output operator, we rewrite Eq. (3) by moving $Lm,\omega $ to the right-hand side to give

where $Rm,\omega =Lm,\omega \u22121$ is the standard resolvent operator. To then specify particular domains for both the response and forcing, we replace $nm,\omega $ with $Bfm,\omega $, where $fm,\omega $ represents an (unrestricted) nonlinear forcing, to the above as

and define the output variable

where ** B** and

**are input and output matrices. The latter matrix,**

*C***C**, is used to isolate the acoustics in the near-field, or propagate fluctuations to the far-field. Each of these cases are detailed in Appendix A. Inserting Eq. (5) into Eq. (6) gives the input-output relationship,

where $Hm,\omega =CRm,\omega B$ is the resolvent input-output operator from $fm,\omega $ to $ym,\omega $. Then by introducing the compressible energy norm of Chu (1965),

(where Ω is the domain volume and superscript $*$ denotes the complex conjugate transpose) *via* the matrix ** W** to the forcing and response ($Wf=Wy=W$) the weighted resolvent input-output operator

is obtained. Resolvent modes may then be found by taking the singular value decomposition of the weighted resolvent input-output operator giving

where the optimal response and forcing modes are contained in the columns of $Um,\omega =Wy\u22121/2U\u0302m,\omega $, with $Um,\omega =[um,\omega 1,um,\omega 2,\u2026,um,\omega N],\u2009Vm,\omega =Wf\u22121/2V\u0302m,\omega ,\u2009Vm,\omega =[vm,\omega 1,vm,\omega 2,\u2026,vm,\omega N]$, and $\Sigma m,\omega =diag(\sigma m,\omega 1,\sigma m,\omega 2,\u2026,\sigma m,\omega N)$ are the optimal gains Towne *et al.* (2018). The unweighted resolvent input-output operator may then be recovered as

### C. Statistics

The statistics we are interested in are contained within the cross-spectral density tensor, which may be found for the desired output space by multiplying the resolvent equation by its complex conjugate transpose and taking the expectation (Towne *et al.*, 2018)

giving

where $Syy,m,\omega $ and $Sff,m,\omega $ are the CSD tensors of the response and the forcing, respectively. For brevity, we drop the subscripts *m* and *ω* and note that all CSD tensors and resolvent matrices must be defined for specific *m* and *ω* pairs in the remainder of the manuscript.

As mentioned in the introduction, this representation shows that if the forcing CSD tensor is known, then the resolvent operator reconstructs the response statistics. However, the forcing CSD is generally unknown. There are at least two potential avenues for modeling it. The first is to directly model $Sff$. To aid in such modeling efforts, $Sff$ may be computed directly from full LES data (Towne *et al.*, 2017), or estimated from limited flow statistics Towne *et al.* (2020). A second approach is to modify the resolvent operator by supplementing the governing linearized equations with an appropriately linearized turbulence model. In Pickering *et al.* (2021), an eddy-viscosity model was considered and LES data were used to determine an optimal eddy viscosity field that would align, insofar as possible, the modes of $Sqq$ (i.e., the full response statistics) with those of $RR*$ (identical to $HH*$ when $C=B=I$). They found this to substantially reduce the magnitude of the off diagonal terms of $Sff$, at least as the near-field coherent structures were concerned, consequently simplifying the number of terms that must be modeled.

In this study, we combine both modeling approaches. We first utilize the eddy-viscosity approximation of Pickering *et al.* (2021) and then estimate a low-order approximation of the forcing CSD for the acoustic field. To do the latter, we return to Eq. (13) and expand the resolvent input-output operator through its singular value decomposition,

and define a covariance matrix $S\beta \beta =V*WfSffWfV$, where $\beta =V*Wff$ is the projection of the forcing upon the resolvent input modes. This gives

which can be rearranged to solve for the covariance matrix,

In its current state, the covariance matrix is exact, maintaining a full size of the system with approximately 10^{11} degrees of freedom (i.e., $S\beta \beta \u2208C5NxNr\xd75NxNr$). To obtain a low-rank model of $S\beta \beta $ from the LES data, we compute $S\beta \beta $ with a truncated set of *n* resolvent modes, $U\u0303\u2208C5NxNr\xd7n$, as,

This reduces the size of the covariance matrix to *n* ×* n*, drastically reducing the number of degrees of freedom to $O(100\u2212101)$.

With $S\u0303\beta \beta $, we may ask several questions: How well does $S\u0303\beta \beta $ reconstruct $Syy$ in the truncated resolvent basis, where the reconstructed CSD is computed as

May $S\u0303\beta \beta $ be further reduced (e.g., neglect off diagonal terms) and can $S\u0303\beta \beta $ be modeled? These questions are addressed in the remainder of the manuscript.

## III. RESULTS

### A. Near acoustic field

We begin by providing detailed results for a single frequency and azimuthal wavenumber pair of the $Mj=1.5$ jet using the RANS eddy-viscosity resolvent operator and a near-field acoustic **C** output matrix [$r/D=[5,20]$ in the fluctuating pressure field, details provided in Eq. (A1)]. Figure 1 presents the first three resolvent modes for $Mj=1.5$, *St* = 0.26, and *m =* 0 computed with the restricted acoustic output domain and then recast in the full domain by setting $C=I$ and calculating $Uq=L\u22121Vy$. The associated gain of these modes, normalized by the first resolvent gain, are $[1,0.17,0.15]$ (and slowly decreasing with higher modes), indicating the first resolvent mode has at least six times the amplification as the following resolvent modes.

The resolvent response modes show a particular pattern of acoustic beams. For the first mode, there is a single, energetic beam, propagating at a shallow angle to the jet axis. The first suboptimal mode consists of two beams, similar to what was found by Jeun *et al.* (2016). This pattern is shown by the next suboptimal mode, with three beams located at the perimeter of the first suboptimal. Although not shown, this behavior continues for further suboptimal modes.

Figure 2 compares three specific (but randomly chosen) realizations of the *m* = 0, *St* = 0.26 field from the LES, ** q**, to the three-mode reconstructions of these fields found by projection. The reconstructions are found by $q\u0303=U\u0303q\alpha \u0303$, where $\alpha \u0303=U\u0303z+*Wzz$ and

**is an acoustic subset of the LES domain ($r/D=[5,6]$ of the pressure field) and $U\u0303z+*$ is the psuedoinverse projecting the LES domain**

*z***and resolvent output domain**

*z***. Further details for using the psuedoinverse to project resolvent modes on other spaces are provided in Appendix B. From Fig. 2, we see that the three resolvent modes are able to accurately reconstruct the different radiation patterns evident in the LES realizations. Clearly, there is constructive and destructive reinforcement amongst the three resolvent modes in order to produce the LES realizations.**

*y*For a more quantitative assessment of the ability of the resolvent modes to reconstruct the acoustic field, we compute and compare the power spectral density (PSD) of the acoustic field, which is located in the diagonal terms of $Syy$. Specifically, we report the difference between the two PSD in decibels,

at $r/D=6$ in Fig. 3. This is again performed for *St =* 0.26, *m =* 0, but is now averaged over all *k =* 153 realizations. In addition to the three resolvent mode set, results are also shown for 5 and 10 mode sets. With just three modes we see that the peak directivity is well captured, with minor improvements (and diminishing returns) in the off-peak directivity with increasing numbers of modes.

We now extend our comparison to Strouhal numbers ranging from 0 to 1 and azimuthal wavenumbers 0–2 and assess the overall ability of the truncated resolvent basis to reconstruct the acoustic field. Figure 4 compares the PSD from the LES to its *n*-rank resolvent-basis reconstructions with *n* = 1, 3, 5, 10, and 20. The rank-1 model results are similar to those of Sinha *et al.* (2014), who used parabolized stability equations and projected onto the first SPOD-mode at each *St* – *m* pair. However, we show here that once additional modes are included, the reconstructions are substantially improved: the 20-mode model shows close agreement with the LES for all frequencies and azimuthal modes, while even the 3-mode model is qualitatively accurate for *m =* 0 and *m =* 1.

To quantify the error between the reconstructed and LES PSD, we compute an error metric,

reflecting the error in acoustic energy. We present this error in Fig. 5 for reconstructions consisting of different numbers of modes. Shown are the errors for two Mach numbers ($Mj=0.9$ and $Mj=1.5$), and three azimuthal wavenumbers, where the filled symbols for the $Mj=1.5$ case provides a reference between the quantitative measure and the qualitative visualization of Fig. 4. Additionally, in order to assess the utility of the eddy viscosity model, we present results both including the eddy viscosity (filled symbols) and neglecting it (hollow symbols).

We find significant improvements in reconstructing the near-acoustic field when including the RANS eddy-viscosity field when compared to results using a constant turbulent Reynolds number of $ReT=3\xd7104$. We note that previous results (Pickering *et al.*, 2021) only considered RANS eddy-viscosity resolvent models with respect to the dominant near-field hydrodynamic SPOD modes. While the rank-1 models for the $Mj=1.5$ jet are similar with and without the eddy viscosity, the remaining reconstructions show a strong and clear advantage to the adopted eddy-viscosity approach. Particularly as sub-optimal modes are added to the basis, the eddy-viscosity model converges rapidly toward the LES, whereas the turbulent-Reynolds-number model shows little improvement. This result is consistent with our previous findings (Pickering *et al.*, 2021), which showed a more profound effect of the eddy viscosity on sub-optimal modes associated with the Orr-mechanism than on modes associated with the Kelvin-Helmholtz mechanism, where the latter are dominant over most of the frequency-wavenumber space being considered here.

Overall, we find the errors related to those with an eddy-viscosity compare favorably when compared to a number of past studies. In particular, the $Mj=1.5$ case at *m =* 0 presents a rank-3 reconstruction with an error of approximately 20%, while previous studies have reported errors >30% when considering a rank-50 reconstruction for only a single frequency (*St* = 0.33) (Jeun *et al.*, 2016). Our approach and metric considers reconstructions over $St=[0,1]$, rather than one individual frequency. Comparing the two Mach numbers, it is apparent that a larger number of modes are required to reconstruct the near acoustic field of the $Mj=0.9$ jet. For example, about ten modes are necessary to obtain a similar quantitative match as compared to just three modes at $Mj=1.5$. This is consistent with multiple past observations where the $Mj=0.9$ jet possesses non-negligible contributions from suboptimal modes that are correlated, or, as described in the time domain, as being linked *via* “jittering” (Cavalieri *et al.*, 2011), thus requiring many modes to reconstruct the acoustic field (Freund and Colonius, 2009; Towne *et al.*, 2015).

### B. Far-field results

We now extend the eddy-viscosity-enhanced resolvent basis to the far-field, and aim to find the modes that are optimal on an arc 100*D* from the nozzle and a range of polar angles from $\varphi =100\xb0$ to $\varphi =160\xb0$ (where $\varphi =180\xb0$ lies on the downstream axis). The domain is depicted in Fig. 6 and the details of the associated output matrix, **C**, are detailed in Appendix A 2.

Figure 7 presents the magnitude of the first three resolvent modes along the arc for both jets at *St =* 0.26 and *m* = 0. The same three-beam structure apparent in Fig. 1 is evident here, with the dominant one-beam mode peaking at $\varphi \u2248150\xb0$. This progression in beam number and location continues in the higher mode numbers not visualized here. Also, plotted in Fig. 7 are the magnitude of modes found *via* spectral proper orthogonal decomposition (SPOD) of the LES data. These modes, which optimally reconstruct the CSD of the far-field arc, are useful to compare to the resolvent modes since a close correspondence between resolvent and SPOD modes indicates that the resolvent mode forcings are mutually uncorrelated (Towne *et al.*, 2018). Indeed, we see a reasonable agreement between the far-field SPOD and resolvent modes. The amplitudes and exact locations vary slightly, but such close agreement suggests that an uncorrelated model may suffice.

Figure 8 shows the near-field signatures of the dominant three far-field modes plotted in Fig. 7 for the $Mj=1.5$ jet. This plot should be directly compared to Fig. 1, which showed the dominant three near-field modes. Outside the jet, the modes are nearly indistinguishable. Within the jet (along the *x* axis), there are differences that can be associated with the larger hydrodynamic wavepacket imprint left in the near-field modes and missing in the far-field ones, consistent with the empirical results of Towne *et al.* (2015).

We now assess how well the computed resolvent modes reconstruct the PSD of the far-field region across $St\u2208[0.1,1]$ and $m=[0,1,2]$ in Fig. 9 using the same PSD error metric as Fig. 5. For both jets at *m* = 0, the rank-1 resolvent reconstruction provides substantial agreement between the LES, at 30% and 40% error for $Mj=0.9$ and 1.5, respectively. The nonzero azimuthal wavenumbers, with the exception of $Mj=1.5$ and *m* = 1, require many additional modes to achieve error levels comparable to the rank-1 *m =* 0 error. This higher-rank behavior is similar to what was observed when reconstructing with near-acoustic-field modes for non-zero azimuthal wavenumbers in Sec. III B.

### C. A simple fit/model

Considering we may reconstruct (i.e., to 30%–40% error) the far-field acoustics at low-rank, we now ask whether we can define a simple forcing model. One approach would be to propose a form of the forcing cross-spectral density tensor, $Sff$, and project this form onto the resolvent input modes to produce a reduced-order matrix $S\u0303\beta \beta $. Despite some clear trends for the dependence of $Sff$ on mean flow quantities (Towne *et al.*, 2017), there does not yet exist a general form for estimating $Sff$. We investigate here an alternative approach of directly estimating $S\u0303\beta \beta $. That is, we focus on modeling the expansion coefficients rather than the forcing itself.

The estimated covariance matrix $S\u0303\beta \beta $ contains the least square reconstruction of the observed data in terms of both the amplitudes and correlations necessary to force each resolvent mode. Where the forcings are uncorrelated, the estimated $S\u0303\beta \beta $ matrix becomes diagonal and only *n* coefficients (albeit at each azimuthal wavenumber and frequency) require modeling. However, even if the forcing is uncorrelated, minor errors or discrepancies in the data, data-processing, computation of resolvent mode, etc., result in a full $S\u0303\beta \beta $ matrix. Further, as rank increases, the statistical uncertainty in the terms becomes greater, reducing our hope for successful modeling. Thus, we explore whether neglecting off diagonal terms is sufficient for a model, but note that this approach provides no guarantees for success; precisely stated, the approximation is not guaranteed to converge as the number of retained modes is increased (Towne *et al.*, 2018).

To limit uncertainty and prevent over-fitting, we assume that the forcing is uncorrelated (i.e.,diagonal) and that the projection of the data with the first resolvent mode,

possesses the lowest uncertainty. These values for the two jets and three azimuthal wavenumbers are shown in Fig. 10. For the $Mj=1.5$ jet, we see that the forcing amplitudes for *m =* 0 and *m =* 1 fall upon lines of constant slope for the most acoustically significant frequency ranges, $St=0.1\u22120.8$. The *m =* 2 data similarly collapse to a line of constant slope; however, the trend is not as clear. Similar observations also hold for the $Mj=0.9$ jet. We also stress that these curves depend on both the data and the resolvent gains, Σ. Including the gains is crucial to collapsing the observed trends.

We now look to fit the data with simple curves of the form,

For the nonzero azimuthal wavenumbers, the data represent the sum of both the clockwise (*m*) and counterclockwise (– *m*) directions about the round jet, such that *a _{m}* represents $a+m+a\u2212m$, or $2a+m$, since $a+m\u2248a\u2212m$, as either rotation about the jet is of equal probability. The exponent,

*b*is unaffected by this symmetry. Figure 10 provides the lines of best fit, where the fits are computed over the region of $St=0.13\u22120.7$ for the $Mj=1.5$ and $St=0.22\u22121$ for the $Mj=0.9$. This choice of low-frequency cutoff is a compromise between keeping the analysis as complete and general as possible, while accounting for the finite computational domain where low-frequency modes are not fully captured within the streamwise extent of the domain, giving the observed plateau in forcing energy. The choice is further justified by noting that the low frequency cutoffs for each case may be related by adjusting the Strouhal number by Mach number, $St0.9=St1.5\xd71.5/0.9$, meaning each range is associated with the same range of acoustic Strouhal number, $Stc\u221e$. The upper bound for the $Mj=0.9$ case extends to

_{m}*St =*1.17, however, we cap the upper bound to

*St*= 1 as done throughout this manuscript.

Table I provides the fitted coefficients for each jet and azimuthal wavenumber. At present, we do not have any physical interpretations of these fits, other than the obvious fact that the power law gives an expected decrease in energy as frequency increases (and thus the length scales of the structures decrease). We suspect we may find similar curves *via* projection of the resolvent forcing modes with the turbulent kinetic energy or other mean-flow quantities, but leave this for future work.

$Mj=0.9$ . | . | . | . | $Mj=1.5$ . | . | . |
---|---|---|---|---|---|---|

Parameter | m = 0 | m = 1 | m = 2 | m = 0 | m = 1 | m = 2 |

a _{m} | $2.65\xd710\u221211$ | $1.38\xd710\u221211$ | $6.14\xd710\u221212$ | $7.10\xd710\u221210$ | $3.89\xd710\u221210$ | $4.66\xd710\u221210$ |

b _{m} | –5.80 | –3.77 | –3.13 | –2.58 | –1.7 | –1.76 |

$Mj=0.9$ . | . | . | . | $Mj=1.5$ . | . | . |
---|---|---|---|---|---|---|

Parameter | m = 0 | m = 1 | m = 2 | m = 0 | m = 1 | m = 2 |

a _{m} | $2.65\xd710\u221211$ | $1.38\xd710\u221211$ | $6.14\xd710\u221212$ | $7.10\xd710\u221210$ | $3.89\xd710\u221210$ | $4.66\xd710\u221210$ |

b _{m} | –5.80 | –3.77 | –3.13 | –2.58 | –1.7 | –1.76 |

To determine how well such curves predict the data, we use the fitted curves to compute,

where $U\u0303$ represents the truncated resolvent basis to rank-*n*. Additionally, as we cannot expect our methods to have accurately captured such large structures in the finite domain used, we use the piece-wise function

where $Stmin=0.22$ and 0.13 for $Mj=0.9$ and $Mj=1.5$, respectively. Figure 9 provides the reconstruction error for the rank-1 model for both jets and three wavenumbers. For the $Mj=1.5$ jet, the rank-1 model yields a close approximation of the reconstructions at 40% for both *m =* 0 and *m =* 1. For $Mj=0.9$, the *m =* 0 error is about 10% worse than the projections, but still at 40% for a rank-1 model. The *m =* 1 case for the $Mj=0.9$ jet and both *m =* 2 cases perform similarly to the perfect reconstructions, but each of these is relatively poor at reducing the error.

With the rank-1 prediction in hand, we conclude by computing the overall sound pressure level (OASPL) found by

Figure 11 presents calculations of the OASPL from both jets using the 100*D* Kirchhoff surface values from the LES and the rank-1 resolvent model considering $m=[0,2]$ and $St=0.13\u22121$. For the $Mj=1.5$ jet, the values are only close at downstream angles and where the jet is the loudest. From $\varphi =140\u2212155$ the model is within 1–3 dB. A striking difference is that the resolvent model peaks at an angle of about 5 degrees higher than the LES data. Interestingly, a similar experiment and simulation in Brès *et al.* (2017) disagreed by the same angle. Although 100*D* data from the experiment is not available, projecting the resolvent modes onto this data would likely result in better alignment, as shifting the LES data by 5 degrees results in a significantly improved estimate (within $0.5\u2009dB$ from $\varphi =135\u2212150$). However, as the ultimate source of the discrepancy is unknown, we avoid making any corrections to the model based on these observations.

We see similar behavior between the Kirchhoff surface (KS) surface and the resolvent model for the transonic case. The rank-1 $m=[0,2]$ resolvent model presents agreement of the peak OASPL to within 2 dB at peak noise angles. We stress that this result for the $Mj=0.9$ jet is rather surprising as many previous studies, although computed in the near-field, found the acoustic field required many modes to agree within 2 dB (Freund and Colonius, 2009; Towne *et al.*, 2015). This shows that the application of both the KS surface to 100*D* and eddy-viscosity model included in our resolvent analysis significantly reduces the rank of the acoustic jet problem. Further, we note that this transonic jet has been extensively verified by experimental data in the near-field and at $\rho =50D$, and, although we extend the results to 100*D*, the peak angles of the KS and the resolvent model are closely aligned when compared to the $Mj=1.5$ case.

## IV. CONCLUSIONS

We formulated resolvent analysis to serve as an acoustic analogy by relating the near-field resolvent forcing to both the near- and far-field acoustic regions. Leveraging the availability of an LES database, we examined resolvent-based reconstructions of the acoustic PSD for turbulent $Mj=0.9$ and $Mj=1.5$ jets. We represented the forcing cross-spectral density matrix with a truncated set of resolvent modes and approximated the amplitudes of the modes with best-fit expansion coefficients of realizations from the LES acoustic field. We found that models consisting of just a single resolvent mode can accurately reconstruct the acoustic field for the first two azimuthal modes for a $Mj=1.5$ jet and the *m* = 0 azimuthal mode for the $Mj=0.9$ jet. To reconstruct higher azimuthal modes, the resolvent basis must be increased to at least five modes (i.e., *M =* 2 and *m =* 1, 2 for $Mj=1.5$ and $Mj=0.9$, respectively). In both jets, the use of an eddy-viscosity model in the resolvent formulation led to clearly superior results compared to a fixed turbulent Reynolds number.

Based on the ability of the rank-1 reconstructions to describe the PSD, we investigated a simple model to collapse the forcing coefficients to one scaling function per azimuthal wavenumber (and Mach number). We found that a power law representation, with only an amplitude and an exponent, suffices to model the coefficient of the optimal resolvent mode. Fortunately, the first resolvent mode contains much of the acoustic energy, and reductions of the gain for this specific mode [related to the Kelvin Helmholtz (KH) mechanism] are likely to provide the greatest reductions in the peak noise of the acoustic field.

The rank-1 $m=[0,2]$ resolvent models estimate the peak noise to within 2 dB for both the $Mj=1.5$ and $Mj=0.9$ jets at peak noise angles. Further, the ability of the resolvent basis to describe much of the acoustic field with only a handful of modes across multiple Mach numbers, a large range of frequencies, and the acoustically dominant azimuthal wavenumbers is promising. This shows that the resolvent framework already contains the appropriate acoustic functions to describe jet noise. In future work, we will seek a fully predictive model by estimating the forcing coefficients from mean flow quantities available from RANS.

## ACKNOWLEDGMENTS

The authors would like to thank André Cavalieri, Oliver Schmidt, and Georgios Rigas for many productive discussions on topics related to this paper. This research was supported by a grant from the Office of Naval Research (Grant Nos. N00014- 16-1-2445 and N00014-20-1–2311) with Dr. Steven Martens as program manager. E.P. was supported by the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program. P.J. acknowledges funding from the Clean Sky 2 Joint Undertaking (JU) under the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 785303. Results reflect only the authors' views and the JU is not responsible for any use that may be made of the information it contains. The LES study was performed at Cascade Technologies, with support from ONR and NAVAIR SBIR project, under the supervision of Dr. John T. Spyropoulos. The main LES calculations were carried out on DoD HPC systems in ERDC DSRC.

### APPENDIX A: OUTPUT MATRICES, C

##### 1. Near-field acoustic output matrix

For analysis of the near-field acoustics, the output matrix ** C** is chosen to only include pressure, $p\u2032=\rho \u2032T\xaf+\rho \xafT\u2032/\gamma Mj2$, in the region

*x*/

*D*= [0, 30],

*r*/

*D*= [5,20]. Ideally, the LES domain would extend from $r/D=[5,20]$ so that the LES could be directly projected onto the resolvent basis; however, the LES database (i.e., the saved data from the LES) only extends to $r/D=6$. Although one could define an output matrix

**that only includes that surface at $r/D=6$, the resolvent modes may still contain hydrodynamic behavior (unless allowed to propagate further from the jet), thus we use the larger domain to ensure the modes are entirely acoustic. Using the larger domain presents a clear loss of orthogonality in the space represented by the LES domain, which is alleviated by truncating the modes to $r/D=[5,6]$ (after computing the resolvent SVD) and implementing a Moore-Penrose inverse such that a least squares fit of the LES in the resolvent basis can be performed. While previous studies have suggested the use of a filter based on the turbulent kinetic energy of the jet within the input matrix**

*C***(Towne**

*B**et al.*, 2017), we take

**to be identity for both the near- and far-field analyses for the sake of generality.**

*B*##### 2. Far-field acoustic output matrix

To define an input-output relationship from the near-field forcing to the far-field acoustics, we introduce a Kirchhoff surface and apply it as a linear operator. We define three radii: *R* as the radial coordinate of the near-field cylindrical surface, *r* as the coordinate pertaining to the far-field cylindrical surface, and *ρ* representing the distance from the nozzle in spherical coordinates (e.g., $\rho /D=100$ for this study). As described in Sec. II B, the input-output problem is defined as

where the output matrix $CR,\rho $ is the total Kirchhoff operator that maps the near-field cylindrical surface, *R*, to the far-field spherical surface, *ρ*. This operator is linearly composed of many Kirchhoff surfaces, $CR,r$, detailed next.

The cylindrical Kirchhoff operator is comprised of several linear operations to ensure accurate results and is defined as,

where $CR$ is a surface selection matrix ($\u2208\mathbb{R}Nsurface\xd75NrNx$), ** N** is an interpolation matrix from a non-uniform grid to a uniform grid with $\Delta x/D=0.025$ ($\u2208\mathbb{R}Nuniform\xd7Nsurface$),

**is a Tukey windowing matrix (using a taper value of 0.75) that extends over the Kirchhoff surface to reduce spectral leakage ($\u2208\mathbb{R}Nuniform\xd7Nuniform$),**

*T***is a padding matrix extending the uniform grid with a total of $2n$ points (**

*P**n*is set to 15) for computing the upstream and downstream wave propagation, as well as ensuring sufficient accuracy in the transform of the initial surface ($\u2208\mathbb{R}2n\xd7Nuniform$),

**is the discrete Fourier transform (DFT) matrix ($\u2208\mathbb{R}2n\xd72n$), and**

*D***contains the derived Hankel functions of the Kirchhoff surface from Freund (2001), with entries along the diagonal for each azimuthal wavenumber, for a specified radial distance,**

*H**r*, from the surface at

*R*($\u2208\mathbb{R}2n\xd72n$).

However, the above operator only supports one specified radial distance from the cylindrical surface at *R*, and a linear combination of $CR,r$ and a proper selection of streamwise points is required to construct a spherical arc. Thus, the linear expression to construct the total Kirchhoff operator is then

where *x _{i}* represents the streamwise location in the 100

*D*arc and

*r*represents the radial extent to which the Kirchhoff surface must propagate from surface

_{i}*R*to the far-field arc

*ρ*for the respective streamwise location. Points are defined along the arc from $\varphi =100\xb0\u2212160\xb0$ with a resolution of $\Delta \varphi =0.5\xb0$.

### APPENDIX B: NON-ORTHOGONAL PROJECTIONS OF RESOLVENT MODES

The statistical relations presented in Sec. II B are valid when ** U** and

**are orthogonal bases in the same space as**

*V***and**

*y***, respectively. However, in the case of the near-field calculations,**

*f***is defined over a larger space than**

*U***and a pseudo inverse must be constructed to find the least-square solution to the above projections. First, we truncate the output modes**

*y***to the output space**

*U**x*/

*D*= [0, 30] and

*r*/

*D*= [5,6] in the pressure field and define the associated output matrix as $Cz$, where

*z*denotes the new restricted space. Applying $Cz$ to both the LES data and resolvent modes gives the ensemble of realizations

**and resolvent modes $Uz$. In addition to reducing the domain space, we also truncate the resolvent response basis to a limited set of**

*z**n*modes, as discussed above, represented as $U\u0303z$. There are now two important consequences of reducing the resolvent domain from $Cy$ to $Cz$. The first is a correction to the gain to the domain $Cz$. Since both output domains share identical input modes we have,

and the gain of the new domain is

where, by definition, $ui,y*Wyui,y=1$. The second is a loss of orthogonality. Fortunately, we may still determine a least squares fit of the data by computing the Moore-Penrose inverse of $Wz1/2U\u0303z,\u2009(Wz1/2U\u0303z)+=(U\u0303z*WzU\u0303z)\u22121U\u0303z*Wz1/2$, and projecting it onto the CSD of *z* to estimate $S\u0303\beta \beta $

This approach is similar to the one taken by Towne *et al.* (2020) for assimilating partially observed flow statistics.