Ethan Pickering

Georgios Rigas

Oliver T Schmidt

Denis Sipp

Tim Colonius

Optimal eddy viscosity for resolvent-based models of coherent structures in turbulent jets

Response modes computed via linear resolvent analysis of a turbulent mean-flow field have been shown to qualitatively capture characteristics of the observed turbulent coherent structures in both wall-bounded and free shear flows. To make such resolvent models predictive, the nonlinear forcing term must be closed. Strategies to do so include imposing self-consistent sets of triadic interactions, proposing various source models, or through turbulence modelling. For the latter, several investigators have proposed using the mean-field eddy viscosity acting linearly on the fluctuation field. In this study, a datadriven approach is taken to quantitatively improve linear resolvent models by deducing an optimal eddy-viscosity field that maximizes the projection of the dominant resolvent mode to the energy-optimal coherent structure educed using spectral proper orthogonal decomposition (SPOD) of data from high-fidelity simulations. We use large-eddy simulation databases for round isothermal jets at subsonic, transonic, and supersonic conditions and show that the optimal eddy viscosity substantially improves the agreement between resolvent and SPOD modes, reaching over 90% agreement at those frequencies where the jet exhibits a low-rank response. We then consider a fixed model for the eddy viscosity and show that with the calibration of a single constant, the results are generally close to the optimal one. In particular, the use of a standard Reynolds-Averaged-Navier-Stokes (RANS) eddy-viscosity resolvent model, with a single coefficient, provides substantial agreement between SPOD and resolvent modes for three turbulent jets and across the most energetic wavenumbers and frequencies.

Introduction

Resolvent analysis (also known as input/output analysis) determines a volumetric distribution of forcing in the frequency domain that gives rise, when acting in a timeinvariant flow, to the most amplified linear response, typically measured in terms of its total kinetic energy. It is an important tool in stability and transition analysis [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF][START_REF] Farrell | Optimal excitation of three-dimensional perturbations in viscous constant shear flow[END_REF][START_REF] Schmid | Stability and transition in shear flows[END_REF][START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF], and has more recently been proposed as a reduced-order model of coherent structures in fully-developed turbulence [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF]Hwang & Cossu 2010b). In the latter context, resolvent analysis can be derived by partitioning the Navier-Stokes equations into terms that are linear and nonlinear with respect to perturbations. Such a rearrangement of the equations is exact, and the equations may be explored without recourse to any further modeling. With varying degrees of formality, similar approaches were proposed in the past [START_REF] Malkus | Outline of a theory of turbulent shear flow[END_REF][START_REF] Michalke | Instability of a compressible circular free jet with consideration of the influence of the jet boundary layer thickness[END_REF][START_REF] Crighton | Stability of slowly diverging jet flow[END_REF][START_REF] Butler | Three-dimensional optimal perturbations in viscous shear flow[END_REF]), but increases in computer power that speed up the singular value decomposition (SVD) of the linear operator using direct LU decomposition (multi-frontal algorithms for sparse systems) have allowed a detailed characterization of the resolvent spectrum in several turbulent, canonical wall-bounded (Hwang & Cossu 2010a,b;[START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF]Sharma & McKeon 2013;[START_REF] Moarref | Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels[END_REF]) and free shear flows [START_REF] Jeun | Input-output analysis of high-speed axisymmetric isothermal jet noise[END_REF][START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF].

At those frequencies where the dominant singular value is significantly larger than the subdominant ones (which we refer to as low-rank behavior), the dominant modes are qualitatively similar to coherent modes extracted from data [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF]. However, when the response is not low rank, a non-trivial structure of the nonlinear forcing terms may lead to discrepancies between resolvent and observed modes. Thus, it is necessary to model the nonlinear forcing to attain resolvent analyses that are quantitatively predictive. Previous studies have considered several approaches for modeling the nonlinear forcing in linear analyses. These include empirical models [START_REF] Bechara | Stochastic approach to noise modeling for free turbulent flows[END_REF][START_REF] Tam | Jet mixing noise from fine-scale turbulence[END_REF][START_REF] Cavalieri | Jittering wave-packet models for subsonic jet noise[END_REF][START_REF] Cavalieri | Coherence decay and its impact on sound radiation by wavepackets[END_REF][START_REF] Towne | A statistical jet-noise model based on the resolvent framework[END_REF], estimation given partial statistics of the response [START_REF] Zare | Colour of turbulence[END_REF][START_REF] Towne | Resolvent-based estimation of space-time flow statistics[END_REF][START_REF] Martini | Resolvent-based optimal estimation of transitional and turbulent flows[END_REF], and/or the use of a turbulent, or eddy, viscosity. An eddy viscosity may be motivated by concepts underlying the triple decomposition [START_REF] Reynolds | Stability of turbulent channel flow, with application to Malkus's theory[END_REF][START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments[END_REF], which identifies the Reynolds stresses as acting on the coherent fluctuations (from both the coherent and incoherent fluctuations), even though the phase average used to define the coherent part of the turbulent-viscosity field is ambiguous in unforced turbulent flows. Many studies have applied eddy-viscosity models in the wall-bounded turbulence literature [START_REF] Del Alamo | Linear energy amplification in turbulent channels[END_REF][START_REF] Cossu | Optimal transient growth and very large-scale structures in turbulent boundary layers[END_REF][START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF]Hwang & Cossu 2010a,b;[START_REF] Hwang | Mesolayer of attached eddies in turbulent channel flow[END_REF][START_REF] Vadarevu | Coherent structures in the linearized impulse response of turbulent channel flow[END_REF][START_REF] Hwang | Attached eddy model revisited using a minimal quasi-linear approximation[END_REF] either through implementation of the [START_REF] Cess | A survey of the literature on heat transfer in turbulent tube flow[END_REF] model or by estimating the eddy-viscosity field via the Reynolds stresses and the mean shear rate of strain. Similarly, global stability analyses have applied eddy-viscosity models to identify and/or control forced or self-sustained resonances in transitional and turbulent flows [START_REF] Crouch | Predicting the onset of flow unsteadiness based on global instability[END_REF][START_REF] Meliga | Sensitivity of 2-D turbulent flow past a D-shaped cylinder using global stability[END_REF][START_REF] Mettot | Quasi-laminar stability and sensitivity analyses for turbulent flows: prediction of low-frequency unsteadiness and passive control[END_REF][START_REF] Sartor | Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile[END_REF]Semeraro et al. 2016a;[START_REF] Tammisola | Coherent structures in a swirl injector at Re = 4800 by nonlinear simulations and linear global modes[END_REF][START_REF] Rukes | An assessment of turbulence models for linear hydrodynamic stability analysis of strongly swirling jets[END_REF][START_REF] Oberleithner | On the impact of swirl on the growth of coherent structures[END_REF]. These studies implemented eddy viscosity on an ad hoc basis, citing improved qualitative agreement or improved integrated energy-densities.

In a more quantitative sense, eddy-viscosity enhanced linear models have also proven useful for assimilating known data to reconstruct observed energy spectra and mean-flow quantities. [START_REF] Moarref | Model-based design of transverse wall oscillations for turbulent drag reduction[END_REF] showed that a data-driven, white-in-time forcing could reproduce the DNS-based turbulent energy spectrum and, similarly, [START_REF] Illingworth | Estimating large-scale structures in wall turbulence using linear models[END_REF] could match DNS energy spectra using time-resolved velocity measurements. More recently, [START_REF] Towne | Resolvent-based estimation of space-time flow statistics[END_REF] showed that incorporating an eddy-viscosity model led to accurate estimates of space-time statistics using partially known data from DNS. Finally, Pickering et al. (2020b) used an eddy-viscosity enhanced resolvent model to reconstruct the large-eddy simulation (LES) acoustic field of transonic and supersonic turbulent jets at a significantly lower rank when compared to their non-eddy-viscosity enhanced computations. Other approaches have implemented eddy-viscosity fields to develop self-consistent models, such as [START_REF] Yim | Self-consistent triple decomposition of the turbulent flow over a backward-facing step under finite amplitude harmonic forcing[END_REF] or [START_REF] Hwang | Attached eddy model revisited using a minimal quasi-linear approximation[END_REF], where the former study coupled a harmonically forced, quasi-linear resolvent analysis with RANS equations, citing eddy viscosity as a necessary link between the coherent and incoherent perturbation dynamics.

Although the utility of eddy-viscosity enhanced linear models for turbulent modeling and control has become increasingly apparent, a quantitative assessment of their effect on turbulent structures is lacking; even more, it is unclear which statistics turbulence models should seek to predict. One appealing target is modes educed by spectral proper orthogonal decomposition (SPOD), as these modes optimally reconstruct the turbulent kinetic energy and represent space-time coherent structures [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. In fact, the SPOD has a theoretical connection with resolvent analysis. [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF] showed that if the resolvent forcing modes, at a given frequency and wavenumber, are mutually uncorrelated, then the resolvent response modes are identical to the SPOD modes. Likewise, discrepancies between the SPOD and resolvent modes imply correlated forcing modes. [START_REF] Morra | On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows[END_REF] applied a similar line of thinking by including an eddy viscosity in their resolvent analysis of turbulent channel flow, showing that the resulting resolvent modes were in greater agreement with the SPOD modes educed from high-fidelity simulation data than resolvent analysis using only molecular viscosity. We extend this approach to turbulent jets, but consider a more general framework. The central question we ask is: how well can the inclusion of an eddy-viscosity model in the resolvent operator approximate the correlations of the forcing cross spectral density tensor? In this approach, an ideal model would render any remaining forcing as uncorrelated, meaning that the resolvent and SPOD modes coincide. We therefore define a data-informed variational problem that seeks an optimal eddy-viscosity field that maximizes the projection of the first SPOD mode on the first resolvent mode. We then show that we can achieve nearly optimal projections using standard eddy-viscosity models, including one directly inferred from a corresponding Reynolds-Average Navier-Stokes (RANS) simulation.

The work presented here is also relevant to a broader debate taking place regarding the interpretation of resolvent analysis. Since we can define the resolvent operator from the full nonlinear equations without introducing approximations or closures, it is attractive to proceed without introducing ad hoc models such as eddy viscosity, since we can still consider the framework exact. With a minor caveat (i.e. while exact, the resolvent decomposition is not necessarily unique as it can depend on the choice of dependent variables used to express the governing equations [START_REF] Karban | Ambiguity in mean-flow-based linear analysis[END_REF]), this implies that the forcing terms are physically interpretable (i.e. measurable) quantities. This perspective is, in our opinion, valuable, and may be pursued alongside efforts (such as the present work) aimed at empirically modeling the forcing. However, there is a subtlety that confounds the separation between "exact" and "modeled" resolvent analyses: namely, it may not be possible to compute, with meaningful accuracy, the exact resolvent modes in high Reynolds number flows, particularly when the mean flow is two-or threedimensional. The fine-scale structure of the modes can require resolutions similar to DNS, and inversion of the resulting linear systems for singular value decomposition can be prohibitive. A survey of resolvent analyses conducted to date on multidimensional base flows shows that a variety of regularizations of the resolvent operator have been used to reduce the computational burden. By regularizations, we mean linear modifications to the operator that, whether through physical or numerical justification, provide results that are free from numerical artifacts or which more closely resemble observed quantities. These include the use of eddy-viscosity models (as discussed at length above), fourthorder numerical filters [START_REF] Jeun | Input-output analysis of high-speed axisymmetric isothermal jet noise[END_REF], effective Reynolds numbers [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF], and linear damping [START_REF] Yeh | Resolvent-analysis-based design of airfoil separation control[END_REF].

From a more general perspective, the present work also has a connection to the building of data-augmented turbulence models [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF]. Here, we specifically target the modeling of unsteady features [START_REF] Wang | Investigations of data-driven closure for subgrid-scale stress in large-eddy simulation[END_REF][START_REF] Maulik | Sub-grid scale model classification and blending through deep learning[END_REF]) and the optimal eddy-viscosity fields found, at each frequency-wavenumber pair, which are analogous to field-inversion steps (also based on variational data-assimilation methods, Foures et al. n cells ∆ta∞/D ∆St subsonic 0.4 4.5 × 10 5 1.117 1.03 15.9 × 10 6 0.2 0.049 transonic 0.9 1.01 × 10 6 1.7 1.15 15.9 × 10 6 0.2 0.022 supersonic 1.5 1.76 × 10 6 3.67 1.45 31 × 10 6 0.1 0.026

Table 1: Parameters, sampling rate, and frequency resolution for the LES.

(2014); [START_REF] Parish | A paradigm for data-driven predictive modeling using field inversion and machine learning[END_REF]) that assist machine learning techniques in generating eddy-viscosity models from mean-flow quantities.

We organize the paper as follows. In § 2 we outline the governing equations, resolvent analysis, and SPOD. In § 3 we discuss the optimization framework developed to match, or align, SPOD and resolvent modes, and the specific eddy-viscosity models examined. § 4 provides the resulting resolvent mode shapes found via the four eddy-viscosity models and § 5 analyzes the associated optimal eddy-viscosity fields. In § 6 we show a favorable impact of the eddy-viscosity models on the subdominant resolvent modes and then conclude the analysis in § 7 by assessing the sensitivity of the RANS eddy-viscosity model. In this final section, we ultimately find a frequency independent RANS eddy-viscosity field that performs well for three turbulent jets (i.e. subsonic, transonic, and supersonic) and their most energetic frequencies (St ∈ [0.05, 1]) and azimuthal wavenumbers (m ∈ N ⊂ [0, 5]).

Methods

The LES database, resolvent analysis, and SPOD were described in [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF] and [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. For brevity, we only recall the main details here.

Large Eddy Simulation database

The flow solver Charles was used to compute the LES databases, including subsonic (Mach 0.4), transonic (Mach 0.9), and supersonic (Mach 1.5) cases; [START_REF] Brès | Unstructured large-eddy simulations of supersonic jets[END_REF] contains the details on the numerical method, meshing, and subgrid-models. Experiments conducted at PPRIME Institute, Poitiers, France were used to validate the Mach 0.4 and 0.9 jets [START_REF] Brès | Importance of the nozzle-exit boundary-layer state in subsonic turbulent jets[END_REF]. Table 1 provides a summary of parameters for the three jets considered. Parameters include the Reynolds number based on diameter Re j = ρ j U j D/µ j (where subscript j specifies the value at the centerline of the jet nozzle exit, ρ is density, µ is viscosity) and the Mach number, M j = U j /a j , where a j is the speed of sound. The simulated M j = 0.4 jet corresponds to the experiments in [START_REF] Cavalieri | Wavepackets in the velocity field of turbulent jets[END_REF]; [START_REF] Jaunet | Two-point coherence of wave packets in turbulent jets[END_REF]; [START_REF] Nogueira | Large-scale, streaky structures in turbulent jets[END_REF] with the same nozzle geometry and similar boundarylayer properties at the nozzle exit. Throughout the manuscript, reported results are nondimensionalized by the mean jet velocity U j , jet diameter D, and dynamic pressure ρ j U 2 j . We report frequencies in Strouhal number, St = f D/U j , where f is the frequency.

Each database comprises 10,000 snapshots separated by ∆ta ∞ /D, where a ∞ is the ambient speed of sound, and is interpolated onto a structured cylindrical grid x, r, θ ∈ [0, 30] × [0, 6] × [0, 2π], where x, r, θ are streamwise, radial, and azimuthal coordinates, respectively. Variables are reported by the vector

q = [ρ, u x , u r , u θ , T ] T , (2.1)
where u x , u r , u θ are the three velocity components, and a standard Reynolds decompo-sition separates the vector into mean, q, and fluctuating, q , components q(x, r, θ, t) = q(x, r) + q (x, r, θ, t).

(2.2)

Resolvent analysis

We start with the nonlinear flow equations of the form

∂q ∂t = F (q), (2.3)
where F is the time-independent compressible Navier-Stokes operator (plus continuity and energy). Substituting equation (2.2) for q and separating terms linear in state perturbations, q , to the left-hand side gives ∂q ∂t -A( q)q = f ( q, q ), (2.4)

where

A( q) = ∂F ∂q ( q) (2.5)
is the linearized flow operator (provided in Appendix B) and f contains the nonlinear terms and any additional external inputs (e.g. environmental noise or perturbations at the boundary).

For the round, statistically-stationary turbulent jets we consider, equation (2.4) is Fourier transformed both temporally and azimuthally to the compact expression

(iωI -A m )q m,ω = f m,ω , (2.6)
where ω = 2πSt is the frequency and m represents the azimuthal wavenumber. We can then rewrite equation (2.6) by defining the resolvent operator, R ω,m = (iωI

-A m ) -1 , q m,ω = R m,ω f m,ω , (2.7)
and introduce the compressible energy norm [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (Part I)[END_REF]) via the matrix W ,

q 1 , q 2 E = q * 1 diag T γ ρM 2 , ρ, ρ, ρ, ρ γ(γ -1) T M 2 q 2 rdrdxdθ = q * 1 W q 2 , (2.8)
to the forcing and response, where W = W f = W q . The resolvent modes under this norm are then found by taking the singular value decomposition of the weighted resolvent operator,

Rm,ω = W 1/2 q R m,ω W -1/2 f = Ũm,ω Σ m,ω Ṽ * m,ω , (2.9)
where the diagonal matrix Σ m,ω contains the ranked gains and the columns of

U m,ω = W -1/2 q Ũm,ω and V m,ω = W -1/2 f
Ṽm,ω contain the response and forcing modes, respectively. These modes are orthonormal in the energy norm, equation (2.8), (2.10) and recover the resolvent operator from equation (2.7) as,

U * m,ω W U m,ω = V * m,ω W V m,ω = I,
R m,ω = U m,ω Σ m,ω V * m,ω W . (2.11)
For the resolvent analysis presented here, just as in [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF], the above equations are discretized in the streamwise and radial directions with fourth-order summation by parts finite differences [START_REF] Mattsson | Summation by parts operators for finite difference approximations of second derivatives[END_REF], while the polar singularity is treated as in [START_REF] Mohseni | Numerical treatment of polar coordinate singularities[END_REF] and non-reflecting boundary conditions are implemented at the domain boundaries.

Spectral Proper Orthogonal Decomposition

SPOD, similar to space-only proper orthogonal decomposition (POD) and originally shown by [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF][START_REF] Lumley | Stochastic tools in turbulence[END_REF], determines an optimal (i.e. in terms of energy) set of orthogonal modes to describe a dataset, but unlike space-only POD, produces modes that express both spatial and temporal correlation in the data. Like dynamic mode decomposition, SPOD modes are computed at unique frequencies. However, through appropriate averaging, SPOD naturally ranks modes by energy and optimally accounts for the statistical variability of turbulent flows [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. Thus, the associated SPOD modes provide the ideal measurement tool to assess modes computed via resolvent analysis.

Decomposing the LES database Q, where Q represents the temporal ensemble of perturbations (q ) found by applying the standard Reynolds decomposition, in the azimuthal and temporal dimensions via the discrete Fourier transform gives the decomposed data matrices, Qm,ω . Multiplying the decomposed matrices, at a particular frequency and azimuthal wavenumber, by their complex conjugate give the cross-spectral density

S m,ω = Qm,ω Q * m,ω , (2.12)
to which we solve the SPOD eigenvalue problem presented by [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF][START_REF] Lumley | Stochastic tools in turbulence[END_REF])

S m,ω WΨ m,ω = Ψ m,ω Λ m,ω .
(2.13)

The SPOD modes form the columns of Ψ m,ω , ranked by the diagonal matrix of eigenvalues Λ m,ω = diag(λ 1 , λ 2 , ..., λ N ). The modes are orthonormal in the norm •, • E , and satisfy Ψ * m,ω W Ψ m,ω = I. As a result, expansion of the cross-spectral density tensor gives,

S m,ω = Ψ m,ω Λ m,ω Ψ * m,ω . (2.14)
In this study, we perform all SPOD computations with a Hamming window and realization sizes of 256 snapshots with 50% overlap, resulting in 78 independent realizations. To avoid ambiguity in referring to computed SPOD and resolvent modes, we use the following notation for the rest of the manuscript. First, all computed mode's subscripts m, ω are dropped, but referenced when necessary in the text. Second, ψ n represents the n-th most energetic SPOD mode, while v n and u n denote the resolvent forcing and response, respectively, that provide the n-th largest linear-amplification gain between v n and u n . Finally, we use the notation ψ 1 : u x when referring to specific components of each mode, as shown here with streamwise velocity.

Using SPOD to inform resolvent analysis

As SPOD provides the optimal description of the second order flow statistics, we wish to use this decomposition to inform our resolvent approach to match such statistics. The connection can be made through multiplication of equation (2.7) by its complex conjugate and then applying the expectation operator to present the relation between the CSD tensors of the forcing and response through the resolvent operator,

S qq = E[qq * ] = E[Rf f * R * ] = RS f f R * .
(2.15)

If q is projected onto the SPOD modes and f is projected onto the input resolvent modes, β = V * W f , where the vector β is the projection coefficients, then we may write when the resolvent and SPOD modes are not identical, which is the case in our study, the forcing coefficients are correlated and this correlation must be modeled.

Ψ ΛΨ * = U ΣS ββ ΣU * , ( 2
Rather than pursuing a direct model of the forcing coefficients, we take an alternative perspective that asks whether a modified resolvent operator, R T , can match one or more of the dominant resolvent and SPOD modes. A trivial solution would be to define the operator by the SPOD expansion, i.e. R T = Ψ , but this operator then corresponds to the (discretization of any) general (non-local) linear operator, rather than a specific partial differential equation (PDE). Instead, a practical model can be obtained by posing a modified PDE of the linearized governing equations with one or more unknown coefficients, and then finding the best choice of coefficients such that the resolvent and SPOD modes are optimally matched. We propose such an approach in the next section by exploiting an eddy-viscosity model, and develop an optimization procedure that fits the parameters to align one, or more, of the most dominant resolvent and SPOD modes.

To the extent that the modified resolvent operator achieves alignment of any one of its output modes with a specific SPOD mode, we may directly interpret the corresponding diagonal entry of S ββ as the forcing amplitude, λ β , required to reproduce the SPOD mode amplitude λ, through the resolvent gain, σ 2 . In other words,

λ n = σ 2
n λ β n for any n where u n = ψ n , (2.17) independent of whether the other modes are aligned (as other modes are orthogonal).

Models considered

We now add an eddy-viscosity model to the linearized governing equations (2.4). We follow the ad hoc model used in (amongst other references) Del [START_REF] Del Alamo | Linear energy amplification in turbulent channels[END_REF] and Hwang & Cossu (2010b), which is typically justified by extending eddy viscosity from its traditional use in modeling the mean Reynolds stresses to modeling the effect of the "background turbulence" on the coherent motion.

The perturbation equations including the eddy viscosity are, with the replacement µ → µ eff = µ j + µ T , identical to the original linearized equations, provided one accounts for the (spatial) variability of µ T (equations provided in Appendix B). There remains an unknown forcing that is the residual between the original forcing and the "coherent" part that is is modeled by the eddy viscosity. Unfortunately, the residual forcing no longer possesses its exact physical interpretation as the nonlinear interactions of resolved modes. However, the advantage is that the resulting response modes can significantly reduce the rank of the problem and lead to a residual forcing CSD that is tractable to model when compared to the forcing CSD of the exactly rearranged equations (Pickering et al. 2020b;[START_REF] Towne | Resolvent-based estimation of space-time flow statistics[END_REF].

In what follows, we refer to the modified linear operator with µ T = 0 as A T and note that the operator depends on the chosen field for µ T , which, upon discretization becomes a vector µ T . Since we assume that µ T is steady and axisymmetric, the operators have a similar temporal/azimuthal Fourier transform that we denote A T m .

We now consider four models for the eddy-viscosity field. The first model directly optimizes the eddy-viscosity field to maximize agreement between the dominant resolvent and SPOD modes. The second model fits an eddy viscosity to the LES mean flow by minimizing the residual in the steady RANS equations. The third model uses an independently computed eddy-viscosity field from a RANS k -model. Finally, we consider a simpler constant eddy-viscosity model based solely upon a turbulent Reynolds number.

For brevity, we refer to the modes computed with the above eddy-viscosity models as EVRA (eddy-viscosity resolvent analysis) modes, while modes termed "baseline" refer to those computed by [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF]. We chose this study as a reference for its extensive comparison of resolvent and SPOD modes across all three turbulent jets and many wavenumbers and frequencies. In the baseline study, they chose an effective Reynolds number of Re T = 3 × 10 4 , a value that is an order of magnitude smaller than the molecular Reynolds number, yet not consistent with the expected magnitude of an eddy viscosity (i.e. Re T << 3 × 10 4 ) . Instead, we regard this intermediate value as a regularization of the resolvent operator. Table 2 summarizes the various models investigated.

For exploratory purposes, we find an eddy-viscosity field that best matches the (so modified) resolvent operator to the measured SPOD modes independently for each frequency and azimuthal mode. The purpose is to gauge the sensitivity of the eddy viscosity value needed to model the different frequencies and azimuthal modes, and should not be interpreted as a proposal for a frequency-dependent eddy viscosity.

Parenthetically, within the following optimization framework we can consider any turbulence model or regularization based on mean-flow quantities. A further example is given in appendix A, where we consider a linear damping model recently proposed for resolvent analysis of unstable base flows [START_REF] Yeh | Resolvent-analysis-based design of airfoil separation control[END_REF]. In the appendix, we find that the linear damping improves agreement, but the performance is generally inferior to the eddy-viscosity models. This is likely due to the monolithic damping effect over all regions and wavenumbers, whereas the eddy-viscosity methods directly addresses the effect of the Reynolds stresses both in regard to specific regions of the flow and the approach's ability to account for spatial gradients in the eddy-viscosity field.

Optimal eddy-viscosity field

Here we develop an optimization, computed independently for each frequency and azimuthal mode, that finds the eddy-viscosity field that is optimal (i.e. the upper bound) in matching the leading resolvent and SPOD modes. To find the analytical expression that determines the sensitivity of mode agreement to an eddy-viscosity field, we use a Lagrangian technique analogous to [START_REF] Brandt | Effect of base-flow variation in noise amplifiers: the flat-plate boundary layer[END_REF] that accounts for the non-modal behavior of the resolvent operator. This technique couples constraints from the governing equations, resolvent analysis, a normalization, and a cost function (agreement of the leading SPOD and resolvent modes), into a Lagrangian functional for whose stationary point provides the desired maximum.

To build the Lagrangian functional, we begin with the forward equation (2.6) and substitute L with L T , the linear operator that includes an eddy-viscosity model. The singular value / singular vector (v 1 , u 1 , σ 1 ) as defined in (2.11) is a solution of both the forward equation (2.6),

v 1 = L T u 1 , (3.1)
where v 1 replaces f as the forcing and u 1 replaces q as the associated response, and the resolvent eigenvalue problem,

W u 1 = σ 2 1 L * T W v 1 . (3.2)
The above resolvent eigenvalue solution is found by taking the energy norm of equation (2.7) and dividing by the forcing energy to give

u * 1 W u 1 v * 1 W v 1 = σ 2 1 = v * 1 R * T W R T v 1 v * 1 W v 1 . (3.3)
Rearranging and eliminating v * n we arrive at

R * T W R T v 1 = σ 2 1 W v 1 , (3.4)
where replacing R T v 1 with u 1 and multiplying both sides by R - * T = L * T recovers equation (3.2). Finally,we define a normalization constraint via,

u 1 , u 1 E = u * 1 W u 1 = 1. (3.5)
The last component of the Lagrangian functional is the cost function,

J = u * 1 W ψ 1 ψ * 1 W u 1 -l 2 µ * T M µ T , (3.6)
where the first term, representing the primary objective, measures the squared projection, what we term the alignment or agreement, between the dominant SPOD mode, ψ 1 , and the first resolvent mode, u 1 . The alignment measure, u * 1 W ψ 1 , is squared to ensure the cost function is real. For brevity, we denote the outer product of the dominant SPOD mode as

Ψ 1 = ψ 1 ψ * 1 = Ψ * 1 .
The cost function may also consider multiple resolvent/SPOD modes by considering a (weighted if desired) sum of the squared alignment terms.

The second term, -l 2 µ * T M µ T , is a Tikhonov regularization that penalizes values of µ T that do not affect the alignment (high values of µ T diminish the value of J ), with M representing the cylindrical quadrature weights of the grid. As done in standard regularization methods, the value of l 2 is chosen high enough to remove the values of µ T in insensitive regions, but also sufficiently small to not interfere with the primary objective [START_REF] Hansen | The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems[END_REF]. This penalization is effective at minimizing the eddy viscosity in non-turbulent regions of the flow such as the far field. A substantial range of l 2 values (i.e. multiple orders of magnitude) removes negligible regions of the eddy-viscosity field from the initial field without an observable drop in the primary objective, alignment between u 1 and ψ 1 .

We now formally construct the Lagrangian functional to include the cost function (3.6), forward equation (3.1), the resolvent eigenvalue problem (3.2), and the normalization constraint (3.5) to give,

L = u * 1 W Ψ 1 W u 1 -l 2 µ * T M µ T -ũ * 1 (L T u 1 -v 1 ) -ṽ * 1 (W u 1 -σ 2 1 L * T W v 1 ) -σ1 (u * 1 W u 1 -1) + c.c., (3.7)
where ( ũ1 , ṽ1 , σ1 ) are Lagrange multipliers, σ1 is real-valued (as the corresponding constraint is real), and c.c. is the complex conjugate. This results in a functional that depends on seven variables,

L([u 1 , v 1 , σ 1 ], [ ũ1 , ṽ1 , σ1 ], µ T ). (3.8)
Algorithm 1 Optimization 1: Initialize. Choose an initial eddy-viscosity/turbulence model and target SPOD mode.

2: while dJ /dµ T = 0 do 3:

Compute the EVRA mode(s).

4:

Solve for the Lagrange multipliers.

5:

Calculate the update direction, dJ /dµ T .

6:

Determine the optimal value of the step α by repeated evaluation of the cost functional along the steepest ascent direction. 7: end while We can find the maximum of the cost function by finding the stationary point of the entire functional (i.e. where variations with respect to each variable are zero). Stationarity with respect to the Lagrange multipliers yields the state equations, which are by definition satisfied, while stationarity with respect to the state variables yields:

∂L ∂u 1 δu 1 = (2W Ψ 1 W u 1 -L * T ũ1 -W ṽ1 -2σ 1 W u 1 ) * δu 1 = 0 (3.9) ∂L ∂v 1 δv 1 = ( ũ1 + σ 2 1 W L T ṽ1 ) * δv 1 = 0 (3.10) ∂L ∂σ 1 δσ 1 = (ṽ * 1 L * T W v 1 ) * δσ 1 = 0, (3.11)
and the condition in the last equation may be simplified into ṽ * 1 L * T W v 1 = ṽ * 1 W u 1 using equation (3.2). The stationary point is subsequently met by constructing the following system of equations and solving for the Lagrange multipliers:

  -L * T -W -2W u 1 W -1 L T σ 2 1 0 0 u * 1 W 0     ũ1 ṽ1 σ1   =   -2W Ψ 1 W u 1 0 0   . (3.12)
The upper left 2 × 2 block is degenerate due to the state equations (3.1) and (3.2) (the couple, ũ1 = W v 1 and ṽ1 = -σ -2 1 u 1 , is in the null-space of this block) and the third column and line regularizes this system. Combining the three equations, one can show that σ1 = u * 1 W Ψ 1 W u 1 , proving that σ1 is a real value. A final variation is taken with respect to the eddy-viscosity, µ T (which may be a scalar or vector quantity), providing the direction of gradient ascent for the eddy-viscosity field,

∂L ∂µ T δµ T = -ũ * 1 ∂L T ∂µ T δµ T u 1 + σ 2 1 ṽ * 1 ∂L * T ∂µ T δµ T W v 1 -2l 2 µ * T M δµ T + c.c (3.13) = dJ dµ T * M δµ T . (3.14)
The gradient at the k th grid point is then:

dJ dµ T k = M -1 km -u * 1,j L * m,ij ũ1,i + σ 2 1 W lj v * 1,l L m,ji ṽ1,i -2l 2 µ T,k + c.c, (3.15)
where

L m,ij = lim →0 L T + δµm,ij -L T ,ij
, δµ m being a null vector except at the m th position where it is equal to 1. This tensor may be obtained either through automatic differentiation of L T with respect to µ T or by finite differences. Full storage of such tensors is not an issue when finite differences, finite volumes, or finite elements are used for the spatial discretization as the resulting tensors are extremely sparse.

The updated optimization parameter is then:

µ (k+1) T = µ (k) T + α dJ dµ T , (3.16)
where k is the iteration number and α is a step size determined through a root finding algorithm or a line search. If multiple SPOD /resolvent modes are considered for the optimization, then one has to solve equation (3.12) for each couple

[Ψ n , (v n , u n , σ n )]
and the total gradient dJ dµ T is the sum of each individual gradient, while the line search for α is performed considering the full cost functional. Although considering multiple modes is theoretically straightforward (and we present one example in § 6), there are two practical issues. Each additional mode brings further complexity to the gradient, increasing computation time, and the quality of SPOD modes, Ψ n , become increasingly noisy with n, thus rendering gains via the optimization as marginal. We discuss the latter issue in more detail throughout the manuscript. Figure 1 presents a schematic of the above optimization framework, including graphical examples from the optimal eddy-viscosity field case at St = 0.6, m = 0, and M j = 0.4.

For some cases, the optimization step imparts a region of negative effective viscosity (i.e. -µ T > µ j ) presenting a challenge in both the physical interpretation and the numerical stability of the resolvent operator. However, negative eddy viscosity is not a unique concept to the algorithm presented. Literature surrounding eddy-viscosity models used in RANS and LES, where the eddy and effective viscosities are identical, attribute physical interpretations of negative eddy-viscosity to backscattering of turbulent energy, which, in many simulations, results in unstable simulations [START_REF] Ghosal | A dynamic localization model for large-eddy simulation of turbulent flows[END_REF]. Common treatment of a negative eddy viscosity has included filtering operations, ensemble averaging in homogeneous directions, and ad hoc clipping of the eddy-viscosity field [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications[END_REF], while inferences of the eddy-viscosity field via a Boussinesq approximation of data are often regularized to remove negative regions (e.g. Semeraro et al. (2016b)). Here, we similarly elect to remove any negative effective viscosity using a simple clipping strategy by setting any negative regions to the molecular value. Although a reduction of effective viscosity below the molecular viscosity is theoretically possible, we found that permitting the optimization to do so either led to numerical instabilities or negligible improvements in alignment.

The topology of the proposed cost function is complex, as µ T involves many degrees of freedom, and our optimizer may return a local rather than global maximum. Therefore, a complete assessment of the sensitivity of initial conditions or demonstration of a global maximum is intractable, but the relative insensitivity of the results to initial guesses and the fact that no other considered method outperforms the full optimization (shown later in figure 4) provides confidence in the robustness of the maxima achieved. For all of the results presented here, we use the optimal constant eddy-viscosity field results (introduced in § 3.4) as the initial condition for the full-field optimizations. Finally, the above optimization is derived considering the full (perturbation) state as the output. The formulation is similar if the input and output spaces are restricted, as shown in Appendix C. Such an extension is of particular use for experiments, or coarse simulations, where observed data may be sparse.

Mean-flow consistent eddy-viscosity model

For many experimental and numerical datasets, including the LES databases used here, an eddy-viscosity field is absent. We circumvent this issue by finding the eddyviscosity field that minimizes the error to which the mean flow satisfies the (zero frequency and axisymetric wavenumber) linearized Navier-Stokes equations, supplemented with an eddy-viscosity model, provided in Appendix B. To do so, we find an eddy-viscosity field that minimizes the residual f given by

L T q = f .
(3.17)

Thus we define the cost function,

J = -f * W f , (3.18)
and develop a Lagrangian functional with the forward equation as the only additional constraint to give

L = -f * W f -ũ * (L T q -f ). (3.19)
Variations with respect to the residual are

∂L ∂f δf = (-2W f + ũ) * δf = 0, (3.20)
and we may directly solve for the Lagrange multipliers as,

ũ = -2W f . (3.21)
Then by taking variations with respect to the eddy-viscosity field gives,

∂L ∂µ T δµ T = -2(W f ) * ∂L T ∂µ T δµ T q. (3.22)
Similar to equation (3.15), we obtain the update step:

dJ dµ T k = -2M -1 km q j L m,ij W il f l , (3.23)
and find the field via a line search. These steps are described in greater detail in the preceding subsection § 3.1. Figure 2 (a) provides the eddy-viscosity field that optimally minimizes the residual of the mean-flow solution. The associated residual field for this model reduced errors to approximately 10% of the original residual field, with the exception where the shear layer is thin near the nozzle. The thin shear-layer region improved by only ≈ 50%, but as shown later in the manuscript, modes in this region are generally less sensitive to the eddy-viscosity field. We refer to this model as the mean-flow consistent eddy-viscosity model, and we optimally tune this field at each frequency by introducing the coefficient, c, µ T = cµ T,M ean . Our interest in the value of c is not to propose a functional of its frequency dependence (or assign to it a physical meaning), but to measure and observe the overall variation and help determine whether a frequency independent coefficient might suffice.

RANS-based eddy-viscosity field

We compute steady-state RANS solutions for each case to assess the applicability of the associated eddy-viscosity field for resolvent analysis. For simplicity, we perform the RANS computations in Fluent. The 2D axisymmetric grid extends 40 diameters in the streamwise directions and 20 diameters in the radial direction with grid spacing mirroring that of the interpolated LES grid scaled to be four times finer, giving 3 × 10 5 grid points. We set the inlet boundary conditions to the base-flow profile from the LES simulations and use the standard 2-equation k -model [START_REF] Launder | The numerical computation of turbulent flows[END_REF] for turbulence modeling. Coefficients used for the model are variants of those suggested by [START_REF] Thies | Computation of turbulent axisymmetric and nonaxisymmetric jet flows using the k-model[END_REF], with turbulent viscosity coefficient C µ = 0.0874, dissipation transport coefficients C 1 = 1.4 and C 2 = 2.02, turbulent Prandtl numbers for kinetic energy σ k = 0.324 and dissipation σ = 0.377, and the turbulent Prandtl number P r T = 0.422. However, the standard κ -model provided in ANSYS does not incorporate the [START_REF] Pope | An explanation of the turbulent round-jet/plane-jet anomaly[END_REF] and [START_REF] Sarkar | The analysis and modelling of dilatational terms in compressible turbulence[END_REF] correction terms used in [START_REF] Thies | Computation of turbulent axisymmetric and nonaxisymmetric jet flows using the k-model[END_REF], requiring a calibration of the mean-flow quantities by introducing a scaling constant a to C µ = 0.0874/a, σ K = 0.324/a, and σ = 0.377/a.

RANS mean-flow quantities closely match those of the LES for each of the three turbulent jets using values for a of 1.2, 1.3, and 1.575, for M j = 0.4, 0.9, and 1.5, respectively. While tuning of the constant a to match LES is not in the spirit of obtaining a universal RANS model, we do so here to give the RANS-generated eddy-viscosity field the best chance at being consistent with the LES results from which the SPOD modes were educed. For a full assessment of the accuracy of RANS predictions for turbulent jets, we refer the reader to [START_REF] Thies | Computation of turbulent axisymmetric and nonaxisymmetric jet flows using the k-model[END_REF]; [START_REF] Georgiadis | Evaluation of modified twoequation turbulence models for jet flow predictions[END_REF].

Figure 2 (b) presents the RANS-predicted eddy-viscosity field for the M j = 0.4 jet, and figure 3 shows near identical agreement with the mean LES streamwise flow. We observe similar agreement in radial velocity, density, and turbulent kinetic energy, and also find close agreement for the M j = 0.9 and 1.5 jets; we do not show these results for brevity. For determination of the optimal RANS-based eddy-viscosity field at each frequency, we take the computed eddy-viscosity fields,

µ T,RAN S = ρC µ k 2 , (3.24)
and introduce the coefficient, c, µ T = cµ T,RAN S (just as in § 3.2). This final relation underscores the difference between the traditional use of eddy viscosity with RANS and ours via resolvent analysis. In the former context, eddy viscosity accounts for all perturbations, while in resolvent analysis, the eddy viscosity is intended to model the effect of nonlinear, triadic interactions and the background turbulence on the linear structures. Thus, a coefficient of c < 1, for resolvent analysis, presents an eddy viscosity that omits a fraction of the overall eddy-viscosity field. As will be shown, we find all optimal values for c to be less than unity. This interpretation may also be applied to the mean-flow consistent eddy-viscosity field presented in the previous sub-section.

Constant eddy-viscosity field

Finally, we consider a simple, constant eddy viscosity, µ T = 1/Re T . We primarily investigate this model because of its use in many turbulent jet studies that used a Reynolds number based either upon the molecular viscosity [START_REF] Jeun | Input-output analysis of high-speed axisymmetric isothermal jet noise[END_REF][START_REF] Lesshafft | Resolvent-based modelling of coherent wavepackets in a turbulent jet[END_REF], on the order of 10 5 -10 6 , or through an effective turbulent viscosity [START_REF] Garnaud | The preferred mode of incompressible jets: linear frequency response analysis[END_REF][START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF], on the order of 10 3 -10 4 . These, quite different, choices inevitably provided discrepancies in amplification gains and mode shapes across each study, particularly at low frequencies (i.e. St < 0.3 for m = 0) -showing that the Reynolds stresses have a substantial impact on resolvent analyses of turbulent jets. Here, we find the optimal Re T at each frequency and azimuthal mode number by a line search.

Optimal SPOD and resolvent mode alignment

In this section, we present modes predicted by the various EVRA models presented in the previous section. We focus on the axisymmetric disturbances, m = 0, for the M j = 0.4 jet, and report results for other azimuthal modes and jet Mach numbers in § 7. We performed optimizations over the frequency range St ∈ [0.05, 1], resulting in the alignment coefficients displayed in figure 4, with alignment defined as |ψ * 1 W u 1 |. This metric not only represents how similar the spatial structures, represented as complex eigenfunctions, are between the dominant resolvent and SPOD modes, but also measures the similarity in distribution of energy amongst the five state variables. A value of 1 signifies perfect agreement, giving both identical agreement in structure and distribution of energy in the state variables. Typically, in this metric, values of approximately 0.4 or greater show qualitative agreement, whereas values less than 0.4 have little visual similarity.

Figure 4 (a) shows that throughout the frequency range considered, the alignments improve considerably from the baseline case (constant eddy viscosity with Re T = 3×10 4 ). The alignment is best for St > 0.3, which corresponds to the frequencies where the jet has a strong, low-rank Kelvin-Helmholtz (KH) response [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF], as highlighted by figure 4 (b), presenting the spectra of the first five SPOD modes and their 95% confidence interval. For this region, St > 0.3, the baseline case gives reasonable (> 75% alignment) results, nonetheless, the eddy-viscosity models still improve the modes to nearly perfect alignment. At lower frequencies, St 0.3, we find the most dramatic increase in alignments, from approximately 10% to 80%. These substantial improvements, at St 0.3, coincide with a change of mode type, from KH to Orr [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF], a viscous, non-modal instability mechanism sensitive to Reynolds number (with rapidly increasing amplification as Reynolds number increases), that dominates the nonoptimized, low-frequency and subdominant regions of the resolvent spectrum for the M j = 0.4 jet. We also find that the optimal eddy-viscosity field provides the greatest alignment among the models, which is at least suggestive that the optimization achieved a global maximum.

Surprisingly, the other eddy-viscosity models produce alignments close to the optimal eddy-viscosity field. The constant eddy-viscosity is nearly optimal at lower frequencies (Orr-type modes), whereas the RANS and optimal mean-flow eddy-viscosity models are more nearly optimal at higher ones. We stress that in the optimal mean-flow, RANS, and constant µ T models, a different optimal value of the coefficient (i.e. c and Re T ) is used at each frequency. We defer a discussion of the sensitivity of these coefficients to § 7.1.

Starting with the lowest frequency, St = 0.05, we now investigate the mode shapes associated with the improved resolvent alignments achieved with the optimized eddyviscosity models. Figure 5 displays the real part of the fluctuating field for all state variables for the dominant SPOD and resolvent modes, comparing resolvent results using both the optimal eddy-viscosity field and the baseline case with constant Re T = 3 × 10 4 . It is immediately apparent that the optimal eddy-viscosity resolvent mode can closely match the observed mode shapes from SPOD for all variables (including the correct distribution of energy), while the baseline resolvent mode bears little resemblance to the SPOD modes for any of the variables.

Despite the increased alignment, there remains an obvious mismatch in u θ between the SPOD and resolvent modes, highlighting a statistical limitation to our approach. For the axisymmetric wavenumber, m = 0, perturbations in the azimuthal velocity must be zero. Both resolvent models meet this constraint, however, the SPOD mode does not. One should then view the nonzero component in the SPOD mode as a statistical error. Compared to the streamwise velocity, u θ is about five times smaller in magnitude, and lacks the coherent wavepacket structure of the other variables. The corresponding u θ contribution in the projection coefficient |ψ * 1 W ψ 1 | is ≈ 0.08, bounding the physical maximum of the optimization to |ψ * 1 W u 1 | 0.92 without considering additional error in the other variables. We link these statistical errors to the weak low-rank behavior with this frequency, where there is little eigenvalue separation between the dominant and subdominant modes [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF]. We may then view the projection-coefficient value of 0.08 as a kind of error bar on the alignments produced by the optimal eddyviscosity field, as it is attempting to align to a mode shape that is (at this frequency) in error by as much as about 10%.

The pressure field, a quantity of particular interest for jet noise, provides a relatively simple representative mode shape for each case. We proceed by visualizing only the fluctuating pressure component for the rest of the study, however, the projection coefficients, |ψ * 1 W u 1 |, account for the full state results. Further, for all response pressure modes presented, we see similar trends and improvements in all flow variables similar to figure 5.

Figure 6 shows the pressure modes at two low frequencies, St = 0.05 and 0.2, and compares the results for all considered eddy-viscosity models. The top row shows the dominant SPOD mode from the LES, the second row gives the dominant resolvent mode for the baseline case, and the remaining rows provide the four optimized models. At low frequencies, the baseline resolvent analysis cannot capture the observed mode shapes, while the optimized eddy-viscosity models have much better alignment with SPOD. The EVRA models increase the projection coefficients by as much as 10-fold and display a wavepacket structure consistent with the SPOD mode. Orr-type modes dominate the low-frequency (i.e. St < 0.3) baseline resolvent spectrum [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF], and we see that the eddy viscosity attenuates these modes in favor of a KH-like response that peaks further upstream, consistent with the observed SPOD modes. Proceeding to higher frequencies, figure 7 displays the dominant fluctuating pressure modes for SPOD and the five EVRA models for St = 0.6 and 1. The baseline projection coefficients are already high for these frequencies, but are further increased with the eddy-viscosity models, reaching 96% for the optimal eddy viscosity. Here the differences in the mode shapes are subtle, with the streamwise extent of the modes shortening from the baseline case to better match the SPOD at both frequencies. At these higher frequencies, the jet response is a clear, low-rank KH wavepacket (a modal, inviscid stability mechanism), and it is thus unsurprising that the results are relatively insensitive to the precise eddy-viscosity model. However, the improved alignment is a product of the non-zero eddy-viscosity field, showing that a turbulence model is still important.

For St = 1, the optimized projection coefficient is falling compared to the St = 0.6 case. This is due to the emergence of Orr-type modes with similar energy as the KH modes. When performing SPOD in limited domains near the nozzle exit, the modal, lowrank KH response continues to dominate at much higher frequencies in the near nozzle region [START_REF] Sasaki | High-frequency wavepackets in turbulent jets[END_REF], but when considering the global response, the KH response becomes inferior, in energy, to the Orr response, which peaks further downstream.

Analysis of the optimized eddy-viscosity fields

The previous section shows that the EVRA approach results in substantial alignment of the dominant resolvent and SPOD modes. In this section, we examine the optimal parameters associated with the eddy-viscosity fields to investigate how the eddy viscosity improved the alignment and to identify potential universalities in modeling coefficients.

Structure of the eddy-viscosity fields

For the constant eddy viscosity, RANS-based, and mean-flow consistent eddy-viscosity fields, the optimization is over a single value, and we plot the optimal values as a function of frequency (still for m = 0) in figure 8 (a,c,d) and the maximum value of the optimal field in 8 (b). We investigated several other metrics for the optimal field and each metric provided similar trends and therefore, we chose µ T ∞ , as it gave the most intuitive comparison against the other scalar quantities. For all models, the frequency dependence of the values is similar, with three regions of interest:

St ∈ [0.05, 0.3], St =∈ [0.3, 0.8],
and St ∈ [0.8, 1].

In the low frequency region, the baseline jet response comprises of spatially extensive Orr-type modes that have a strong Reynolds number dependence, requiring a relatively larger eddy viscosity to damp them. For St = 0.05 the ratio of the optimal effective Reynolds number to the molecular Reynolds number is µ T /µ j ≈ 13, 500, a four orderof-magnitude difference when compared to the molecular viscosity.

In the moderate frequency regime, where the baseline spectrum transitions from the broadband, viscous Orr mechanism to the low-rank, inviscid KH mechanism, eddy viscosity becomes less important, and we expect (confirming below, in § 7.1) insensitivity to the overall value based on the relatively favorable alignment achieved in the baseline case. As frequency increases, the responses transition back to a mix of KH and Orr-type waves, with a progression towards broadband, viscous Orr modes at higher frequency.

At these higher frequencies, we see that the low-frequency dependence on inverse effective Reynolds number resumes, similar to low the frequencies. Interestingly, this trend shows that at higher frequencies Re T → Re j such that the effect of eddy viscosity "turns-off" as frequency increases and the associated wavepacket wavelength becomes small (i.e. approaching finer-scale turbulence), as expected on physical grounds. Similar trends are observed for the mean-flow consistent and RANS eddy-viscosity coefficients. For both eddy-viscosity fields, c is less than unity and only varies from 0.7 to 0.1 with a few exceptions. These values suggest that the optimal eddy-viscosity model is a fraction of the total RANS or mean-flow consistent eddy-viscosity models that integrate all perturbations.

For the full-field eddy-viscosity optimization, we stress that its primary purpose is to determine what may be an upper bound for how well any eddy-viscosity model could perform. Given that the alignments between the resolvent and SPOD modes were not significantly higher for the optimized scheme than for the modeled eddyviscosity approaches (with optimal parameters), the detailed eddy-viscosity fields are of lesser importance. Still, some aspects of the physics, such as the spatial locations where Reynolds stresses become important for each frequency, are apparent in the optimized fields. Figure 9 presents the optimized fields for two selected Strouhal numbers, comparing them to both the RANS and mean-flow consistent eddy-viscosity fields scaled by their optimal coefficient c at each frequency. In addition, the dominant resolvent mode, computed with the displayed optimal-eddy-viscosity field, is shown for comparison with the eddy-viscosity fields. The contours for the eddy-viscosity fields are set from 0 to the maximum value of the St = 0.6 optimal eddy-viscosity field.

Overall, both frequencies present optimal eddy-viscosity fields that are complex, unsurprising given the ability of the optimization to choose any eddy-viscosity field, constrained only by the structure of the equations and positivity. The optimal eddy-viscosity fields pinpoint the locations where linear structures break down (i.e. where nonlinearities/Reynolds stresses become important) and inform what features an eddy-viscosity model must include. In both cases, the optimization removes viscosity from the potential core (i.e. the interior region of the jet relative to the critical layer), when compared to the initial guess, while increasing the turbulent viscosity just outside of the critical layer. The increase in eddy-viscosity is most often observed just downstream of the peak amplitude of the wavepacket, coinciding with each wavepacket's decay downstream.

Although not entirely clear from figure 9, these findings are reasonably consistent with each of the modeled eddy-viscosity fields when restricting the view to the region where the resolvent/SPOD mode has significant amplitude. We can see that both the RANS and mean-flow consistent eddy-viscosity fields present similar features as the optimal field, explaining the ability of each model to achieve nearly optimal results. We will show in the following section how such features also explain the ability of the RANS and mean-flow consistent models to predict the subdominant modes, which require further turbulence modeling downstream.

Alignment of subdominant modes

Although the optimization presented only aligns the dominant SPOD and resolvent modes, subdominant modes are also of interest, particularly as they are necessary to reconstruct flow statistics in the near field and are relevant for modeling coherence decay associated with the "jittering of wavepackets" to produce sound [START_REF] Cavalieri | Jittering wave-packet models for subsonic jet noise[END_REF]. In this section we seek to answer two questions: whether alignment with only the dominant mode substantially alters the alignment of the subdominant modes and the effect of expanding the optimization to subdominant modes. We first assess the former case using the optimal parameters for each method. We show the computed subdominant modes in figure 10 for modes 2 and 3 for the St = 0.6, m = 0 frequency-wavenumber pair.

Comparing the second mode to the baseline case (Re T = 3 × 10 4 ), we find that all EVRA models give significantly improved alignments, reaching ≈70% for the RANS and mean-flow consistent models. Both the RANS and mean-flow consistent models are superior to the optimal eddy-viscosity field, which is only fitted to align the dominant mode. The RANS and mean-field models are also superior for the third, fourth, and fifth modes (the latter two not shown for brevity), but with an alignment that falls off with increasing mode number.

To observe how well the optimization of the first SPOD mode models the forcing statistics (i.e. diagonalizes the forcing CSD S f f ), we compare projections of the first five SPOD modes with the first five modes from each eddy-viscosity method (including a 2mode optimization described next) in figure 11. The plots show that the EVRA models, in particular the RANS and mean-flow consistent models, are superior at diagonalizing the CSD when compared to the baseline case.

Although the optimal eddy-viscosity field, aligned only with the dominant SPOD mode, shows improvements in the subdominant modes, we can extend the optimization to align an arbitrary number of subdominant modes and achieve alignment superior to any eddyviscosity model. However, convergence issues with increasing SPOD mode number suggest that optimizing for many modes (e.g. n > 5) would have marginal returns. For this study, we present only the optimization of both the first and second modes at St = 0.6, m = 0 to show the generality of the optimization framework and the physical implications of the associated eddy-viscosity field for the subdominant modes.

Figure 12 presents the aligned resolvent mode via the optimization and the associated eddy-viscosity field for the first subdominant mode. By including the second SPOD mode, We find that this additional downstream eddy-viscosity, present in both the RANS and mean-flow consistent models, is responsible for the increased subdominant mode alignment. Considering the RANS (and mean-flow) model also shows similar downstream structure, we investigate its merit for a predictive model in the next section.

Towards a predictive EVRA model for turbulent jets

Through the previous sections, we have shown that both the RANS and mean-flow eddy-viscosity models perform well across Strouhal numbers from 0.05 to 1 at m = 0, provided the overall constant associated with their application to the fields is optimal (at each frequency and azimuthal mode number). In this section, we consider the sensitivity of the results regarding the of a (and wavenumber) independent and that over a range of frequencies and azimuthal mode numbers, alignments are relatively insensitive to the choice of a constant, such that a single, universal value may be acceptable. While both RANS and mean-flow consistent models both performed well with optimal coefficients, we focus only on the RANS k -model, as it is better regarded as universal across a range of flows. We then apply EVRA-RANS to the M j = 0.4 jet using a single constant to six azimuthal wavenumbers, m = 0 -5, and find substantially improved predictions when compared to the baseline. We also find similar observations when using the same EVRA-RANS model for both the transonic and supersonic jets. Finally, we present the effect of the eddy viscosity on the resolvent spectra.

Frequency and azimuthal mode sensitivity

The optimal RANS coefficients (figure 8) ranged from c = 0.7 -0.004, with a relatively constant region, c = 0.5, for moderate frequencies and, considering the fully optimized eddy-viscosity field produced only marginally improved alignments for most cases, the results may not be sensitive to the precise constant. We test this hypothesis for the RANS model across a range of frequencies with proposed "universal" values of constant c = [1, 0.5, 0.32, 0.2, 0.08]. We plot the resulting alignments versus frequency in figure 13. With little compromise, compared to the optimal constant for each frequency, a single constant of c = 0.2 provides significant alignment across all frequencies up to St = 1.

Although not shown for brevity, we found similar observations using c = 0.08 -1 for all three Mach numbers and six azimuthal wavenumbers. In these cases, not only did c = 0.2 give the best overall alignment, but the alignments were comparably insensitive to the value of c chosen over this range.

We do not present a rigorous justification for the value of c = 0.2, however, c may have a connection with LES eddy-viscosity modeling. If we compare the empirical coefficients used in this RANS eddy-viscosity analysis, cC µ , to the square of the Smagorinsky coefficient, C 2 s , used for SGS eddy-viscosity fields in LES [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment[END_REF], we can find that the value of c = 0.2 lies within the bounds of wall-bounded and isotropic turbulence. Setting equal the products of each set of coefficients, we find c = C 2 s /C µ . Then using the M j = 0.4 RANS coefficient, C µ = 0.073, and the commonly used range of C s , 0.1 for wall-bounded flows and 0.18 for isotropic turbulence [START_REF] Zhiyin | Large-eddy simulation: Past, present and the future[END_REF], we find that c = [0.14 -0.45] (ranges for the M j = 0.9, 1.5 jet are c = [0.15 -0.48], [0.17 -0.57], respectively). This range approximately corresponds to the acceptable values found in figure 13.

For nonzero azimuthal modes, figure 14 presents the alignment of the EVRA-RANS model with SPOD using c = 0.2 and the baseline case for m = 0 -5. The EVRA-RANS model substantially increases the alignments for all nonzero wavenumbers. The results for m = 1 are particularly encouraging, with a uniform, 80% alignment across all frequencies. Azimuthal modes greater than 1 result in poorer alignment, albeit much improved compared to the baseline case, especially when m > 2.

Expanding to nonzero azimuthal wavenumbers, the eddy-viscosity field also affects a third mechanism observed in the global SPOD spectrum (as St → 0), the lift-up mechanism (Pickering et al. 2020a). Similar to the Orr mechanism, the lift-up mechanism They also observed that, in turbulent jets, the spatial extent of resolvent modes increase as frequency decreases and that without an eddy-viscosity, modes extend indefinitely downstream for St = 0. This is analogous to theory surrounding streaks where the lift-up mechanism presents a rapid spatial growth of streamwise streaks until viscous dissipation becomes dominant and the structures decay [START_REF] Hultgren | Algebraic growth of disturbances in a laminar boundary layer[END_REF]. Considering the significant improvements between alignments for low-frequency and nonzero wavenumbers, we find the lift-up mechanism to also be sensitive to an eddyviscosity model.

Transonic and supersonic turbulent jets

We now generalize the RANS-EVRA model performance for both M j = 0.9 and 1.5 turbulent jets using c = 0.2. Figure 15 provides the alignments across frequencies and azimuthal wavenumbers for each. The transonic jet gives substantial agreement for m = 0 and m = 1 at about 80% for much of the frequency range, while m = 2 gives alignments of 60%, on average. For the supersonic jet, the agreement is not as favorable, however, much improved from the Re T = 3 × 10 4 alignments (not shown here for brevity).

The RANS eddy-viscosity model increases many of the alignments, however, poor alignments remain, and these alignments appear to correspond to SPOD spectra without large energy separation. As shown earlier in figure 4, EVRA and SPOD modes aligned best when there exists large eigenvalue separation between the first and second SPOD mode. We find similar behavior here for all cases. Figure 16 presents the SPOD spectra of the first five modes across all six azimuthal wavenumbers and three turbulent jets, with their associated 95 % confidence intervals in light blue. A handful of the spectra show a clear separation between mode energies, such as those between the first and second mode for M j = 0.4, m = 0 and 1, for M j = 0.9, m = 0 and 1 (and higher frequencies for m = 2 -5), and for M j = 1.5, m = 1. In each case where there is large eigenvalue separation, we find, from figures 14 and 15, significantly greater agreement in projection coefficients between the resolvent and SPOD modes, while finding poor projections for cases without clear separation in eigenvalues. We also observe this for the subdominant modes investigated in the M j = 0.4, m = 0 case in § 6. These observations point to a limitation to our method when comparing EVRA modes with SPOD modes. For the SPOD modes without clear eigenvalue separation, the eigenvalues themselves fall within the uncertainty bands (i.e. 95 % confidence interval) of the other modes. The eigenvectors corresponding to these eigenvalues are expected to have, at best, similar uncertainty levels. Thus, without more data, it is not possible to attribute the lack of agreement to a failure of the EVRA ansatz.

Singular values

We return to the m = 0, St ∈ [0.05, 1] case to assess the EVRA-RANS c = 0.2 model's effect on the singular values and compare them to the baseline case and the SPOD eigenvalues. Figure 17 provides the spectra of the first five modes for SPOD (accompanied by a shaded region providing the 95% confidence interval of the eigenvalues), the baseline resolvent model, and the RANS-EVRA model (using c = 0.2) for m = 0. Comparing the resolvent spectra to the SPOD spectra, we immediately see that the separation between λ (i.e. the ratio between λ n /λ n+1 ) and σ 2 of either resolvent models does not compare favorably. In fact, the RANS-EVRA spectra has increased its energetic separation when compared to the baseline case.

This behavior may be linked to multiple (in this case two for m = 0) distinct mechanisms represented in the flow, the KH and Orr-mechanisms. As detailed earlier, the inclusion of an eddy-viscosity model presents a substantial effect on the Orr modes significantly reducing the streamwise extent of each mode, while the KH modes are relatively unchanged. We observe an analogous effect here in figure 17 where the singular values related to the Orr mechanism decrease substantially, pulling away from the unaffected singular values of the KH mechanism, resulting in a much larger separation between singular values than is observed between the SPOD eigenvalues. This sensitivity of Orr modes to an eddy-viscosity was also observed in [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF] at St = 0.6, m = 0 when adjusting Re T , finding that the squared singular values of the subdominant Orr modes scaled as Re 1.2 T . We observe the same effect using the RANS eddy-viscosity model, interestingly (and perhaps unsurprising given the preceding discussions), figure 17 (c) provides similar values as those reported by [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF] at St = 0.6, m = 0 when using Re T = 10 3 .

Figure 17 (a) and (c) also show that the forcing amplitudes, λ β , are not uniform in turbulent jets, contrary to a customary assumption used in resolvent analysis where Λ β = αI, with α as an arbitrary constant [START_REF] Morra | On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows[END_REF][START_REF] Lesshafft | Resolvent-based modelling of coherent wavepackets in a turbulent jet[END_REF][START_REF] Hwang | Attached eddy model revisited using a minimal quasi-linear approximation[END_REF]. Focusing on only the first and second resolvent and SPOD modes for St = 0.6 and m = 0, where mode alignments are 95% and 69%, respectively (the optimal-field case increases the latter value to 77% without appreciably changing the singular value), we may assume that the diagonal components of the forcing, λ β,1 and λ β,2 , account for nearly all the energetic contributions by these two modes. As shown by equation (2.17), this assumption allows for a one-to-one comparison between the first two SPOD eigenvalues, resolvent singular values, and forcing amplitudes (i.e. λ β,n = λ n σ -2 n for n = 1, 2). If the customary assumption of uniform forcing is applied, Λ β = αI, then α = λ 1 σ -2 1 = λ 2 σ -2 2 or, alternatively, λ 1 /λ 2 = σ 2 1 /σ 2 2 , and figure 17 shows this cannot be true. Therefore, unless the SPOD and resolvent spectra are equivalent, we must model or estimate the non-trivial forcing amplitudes.

The sizeable difference between the singular values reflects the forcing of different mechanisms at significantly different amplitudes in the flow. Pickering et al. (2020a) showed that there are three distinct spatial regions that lead to the most efficient amplification of the KH, Orr, and the lift-up mechanisms. They found that regions localized near the nozzle where perturbations are smaller, associate with KH-type responses, while regions downstream and near the end of the potential core where perturbations are significantly larger, support Orr-type responses. Considering these observations, a logical next step in completing a resolvent-based turbulence model is to tie the forcing amplitude of different modes to the turbulence intensities in the respective regions that force them.

Conclusions

We developed a data-informed optimization that quantitatively tested the extent to which an eddy-viscosity model improves the alignment (i.e. agreement) between observed large-scale structures, educed via SPOD, and those computed from resolvent analysis. This eddy-viscosity approach acts as a proxy for modeling the effect of turbulence on large-scale structures, and we found this approach provides substantial improvements in agreement (i.e. when compared to a baseline case that used a constant eddy-viscosity model corresponding to a value of Re T = 3 × 10 4 ). By directly optimizing the eddy viscosity field to achieve the best alignment, we found alignments between resolvent and SPOD modes as high as 96% or improvements of over ten-fold from the baseline alignment (i.e. 8% to 80%).

Across the frequencies and wavenumbers considered, the addition of an eddy-viscosity model to the resolvent operator highlighted its effect on the different amplifications mechanisms in the turbulent jet: Orr-type, KH-type, and lift-up. Although eddy-viscosity models improved modes related to the KH-type mode, we found KH modes to be rather insensitive to the eddy-viscosity field, a result expected from the inviscid nature of the inflectional KH instability. For resolvent modes associated with the Orr and lift-up mechanisms, known to arise from nonlinear interactions, we found significant sensitivity. Resolvent modes computed without a sufficient eddy-viscosity model were visually unrecognizable from their SPOD counterpart, while those computed with an eddy-viscosity model aligned to nearly 80%.

The optimal eddy-viscosity field also provided an upper bound for mode agreement, providing a benchmark to assess three additional eddy-viscosity models. Of these models, we found that traditional eddy-viscosity models (e.g. RANS based) perform nearly as well as the optimal eddy-viscosity models in aligning the most energetic mode. The traditional models even outperformed the optimal model (i.e. optimal in the first mode) when considering the subdominant modes, giving the greatest diagonalization of the forcing CSD at m = 0, St = 0.6 (i.e. ability to model the effect of nonlinear forcing), leading to a more efficient resolvent basis for describing turbulent jets.

Finally, we tested the modeling potential of a RANS-inferred EVRA through a sensitivity analysis and observed its performance over frequency, azimuthal wavenumber, and Mach number. We found the sensitivity of the RANS-based EVRA model's calibration constant, c, to be weak, giving similar agreement for coefficients ranging over an order of magnitude. We found that the coefficients resulting in similar agreement fell within a range of values that may have a connection to the Smagorinsky coefficient for sub-grid scale modeling of eddy viscosity in LES. Choosing a frequency-independent RANS-EVRA model (i.e. c = 0.2), we tested its performance across six azimuthal frequencies and three turbulent jets, spanning subsonic, transonic, and supersonic regimes. For the first three azimuthal wavenumbers (i.e. m = 0 -2), we observed substantially increased alignments for all three turbulent jets and across Strouhal numbers St ∈ [0.05, 1]. Overall, these results show that "classical" eddy-viscosity models (RANS or a mean-flow consistent model) aid in estimating the impact of the Reynolds stresses for resolvent analysis.

While the present data-driven analysis points to the efficacy of relatively simple eddyviscosity-based models for modeling the effect of nonlinear forcing and providing a more efficient resolvent basis, there remains a need for refinements to this approach and careful comparison and consideration of alternative formulations. whose solution provides the Lagrange multipliers, ũ1 , ṽ1 , σ1 . A difficulty that arises in building equation C 5 is that the term L -1 T is a large, dense matrix. When L, B, C, and the weighting matrices are sparse, we may instead introduce auxiliary variables through

L T η1 = B ũ1 L *
T ζ1 = C * ṽ1 , whereupon equation C 5 may be written as a larger, but now sparse, system of equations

      L * T 0 -C * 0 0 0 L T 0 -B 0 0 W * y C I 0 W y u 1 B * 0 0 σ 2 1 W * f 0 0 0 0 v * 1 W f 0             ζ1 η1 ũ1 ṽ1 σ1       =       0 0 2W y Ψ 1 W y u 1 0 0       . (C 6)
The above presents a general optimization framework for aligning/matching any inputoutput resolvent analysis to data (i.e. here we use SPOD modes, but Ψ need not be restricted to SPOD modes). Variations with respect to any parameter of the resolvent operator may now be made to investigate their effect on modelling (or assimilating) known quantities.

Figure 1 :

 1 Figure1: Schematic of the optimization framework for determining the optimal eddyviscosity field that maximizes the alignment/agreement between computed resolvent modes, u 1 , and educed SPOD modes, ψ 1 . Included graphics are from implementation of the full-field eddy-viscosity model at St = 0.6, m = 0, and M j = 0.4.

Figure 2 :

 2 Figure 2: (a) Mean-flow consistent eddy-viscosity model computed at zero frequency and azimuthal wavenumber. (b) Eddy-viscosity field computed via a RANS simulation for the M j = 0.4 jet, c = 1.

Figure 3 :

 3 Figure 3: Mean-flow profiles of both the M j = 0.4 LES and RANS, where the RANS simulation was tuned to best match the LES mean flow. (a) presents the streamwise mean velocity at three radial locations, r/D = 0.25, 0.5, 1, versus streamwise distance from the nozzle, while (b) gives the streamwise mean velocity at three streamwise locations, x/D = 0.5, 5, 10, versus radial distance.

Figure 4 :

 4 Figure 4: (a) Optimal alignments for all methods investigated including the baseline case, Re T = 3 × 10 4 . (b) SPOD eigenvalue spectra of the first five modes for m = 0, including the 95% confidence intervals and the modes associated with the Kelvin-Helmholtz and Orr mechanisms.

Figure 5 :

 5 Figure 5: Real component of the fluctuating response state variables, q = [ρ, u x , u r , u θ , T ], and pressure, p, at St = 0.05, m = 0. The columns display SPOD (ψ 1 ), optimal eddy viscosity (u 1 ), and baseline (u 1 ) modes from left to right, respectively. Contours ( ) are given by ±0.5||ψ 1 : •|| ∞ of the SPOD mode, where • is the fluctuating variable in question (with ||ψ 1 : •|| ∞ values: [ρ, u x , u r , u θ , T, p] =[2.8, 198.6, 46.0, 37.2, 1.2, 10.4] × 10 -3 ).

Figure 6 :

 6 Figure 6: Real component of the response pressure fluctuations ( , ±0.5||ψ 1 : p|| ∞ ) for St = 0.05 and St = 0.2 in the left and right columns, respectively. Row 1 presents the dominant SPOD mode for which the optimization seeks to match. The following rows present results for the baseline, optimal eddy-viscosity field, mean-flow consistent model, RANS eddy-viscosity model, and the optimal turbulent Reynolds number.

Figure 7 :

 7 Figure 7: Real component of the response pressure fluctuations for St = 0.6 and St = 1 in the left and right columns, respectively. Rows present the equivalent methods as described in figure 6.

Figure 8 :

 8 Figure 8: The optimal parameters across St ∈ [0.05, 1] for (a) the optimal constant field 1/Re T , (b) optimal eddy-viscosity field model, (c) the mean-flow consistent model, and (d) the optimal RANS model. The optimal eddy-viscosity field parameter shown is the maximum value of the field at each frequency, µ T ∞ , while the latter two models present the optimal coefficient c. The associated alignments for each model/parameter are shown in figure 4.

Figure 9 :

 9 Figure9: Comparisons of the optimal eddy-viscosity fields (i.e. full-field optimal, meanflow consistent, and RANS) and the associated dominant resolvent mode found via the optimization for St = 0.2 and 0.6. Contours for all six eddy-viscosity fields are set from 0 to 3 × 10 -3 .

Figure 10 :Figure 11 :

 1011 Figure 10: Subdominant modes 2 and 3 at St = 0.6, m = 0 in the left and right columns respectively for SPOD, baseline, and all EVRA models.

Figure 12 :Figure 13 :

 1213 Figure12: The second subdominant mode at St = 0.6 and the associated eddy-viscosity field that provides the optimal alignment for both modes. The contour for the eddyviscosity field is set to the same value as those shown in figure9from 0 to 3 × 10 -3 .

Figure 14 :

 14 Figure 14: Alignments for frequencies, St ∈ [0.05, 1], and azimuthal wavenumbers, m = 0 -5, for the (a) RANS eddy-viscosity model using c = 0.2 and the (b) baseline, constant eddy-viscosity case (i.e. Re T = 3 × 10 4 ).

Figure 15 :

 15 Figure 15: Alignments using the RANS eddy-viscosity model with coefficient c = 0.2 across Strouhal numbers St ∈ [0.05, 1] and azimuthal wavenumbers m = 0 -5 for the (a) M j = 0.9 and (b) 1.5 jets.

Figure 16 :

 16 Figure 16: Spectra, and their associated 95 % confidence interval in light blue, of the first five SPOD modes for azimuthal wavenumbers m = 0 -5 from left to right and the subsonic, transonic, and supersonic jets from top to bottom, respectively.

Figure 17 :

 17 Figure 17: Spectra of first five (a) SPOD, (b) baseline resolvent, and (c) the RANS eddy-viscosity model resolvent modes at m = 0 for St ∈ [0.05, 1].

Figure 18 :

 18 Figure 18: Optimal alignments for the linear damping term and the baseline case, Re T = 3 × 10 4 ..

  .16) which highlights that if the forcing coefficients are uncorrelated (S ββ = Λ β ) then the resolvent modes would be equivalent to the SPOD modes[START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. Conversely,

	Turbulence model	µT form	Optimal parameter LES data used Abbreviation
	Baseline*	1/ReT = 3. 3×10 -5	-	-	Baseline
	Optimal field	µT (x)	µT (x)	Ψ	Opt. µT
	Mean-flow consistent	cµT (x)	c	q	Mean µT
	RANS	cρCµk 2 /	c	-	RANS µT
	Turbulent Re	1/ReT	1/ReT	-	ReT,Opt

Table 2 :

 2 Turbulence models investigated in this study. The baseline* case refers to the results of[START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF].

Acknowledgments

This research was supported by a grant from the Office of Naval Research (grant No. N00014-16-1-2445) 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. 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.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Linear damping term

Besides the studied eddy-viscosity models, we also investigated the impact of a linear damping term, which is equivalent to a finite-time-horizon resolvent analysis introduced by [START_REF] Jovanović | Modeling, and control of spatially distributed systems[END_REF], recently studied by [START_REF] Yeh | Resolvent-analysis-based design of airfoil separation control[END_REF] to localize the resolvent forcing and response modes on an airfoil. For this model, we modify the operator so that,

where β = 1/τ > 0, and τ is the desired temporal decay rate. We then find the value of β that best matches the dominant resolvent and SPOD modes.

Figure 18 presents the agreements/alignments for the linear damping case. Although linear damping improves alignments, the performance is significantly inferior to the eddyviscosity models, likely because of its monolithic damping effect over all wavenumbers, whereas the eddy-viscosity methods directly address the effect of the Reynolds stresses. Considering its suboptimal performance when compared to eddy-viscosity models, we only present results for the M j = 0.4, m = 0, and St ∈ [0.05, 1] cases.

Appendix B. Governing equations

Conservation of mass, momentum, and energy for a compressible, Newtonian fluid are written as,

respectively, where Θ = ∇ • u is the dilatation. We take P r ∞ = 0.7 and γ = 1.4 as constants. The equations have been made nondimensional with the jet density (ρ j ), speed (U j ) , and diameter, D. The nondimensional viscosity, µ = 1 Rej , is also a constant. Applying the Reynolds decomposition (i.e. q(x, t) = q(x) + q (x, t)) to the above equations and separating terms that are linear and nonlinear in the fluctuations to the left-and right-hand sides, respectively, gives

, and where we have grouped all the nonlinear terms as forcing terms on the right-hand sides.

The left-hand side is then transformed to a cylindrical coordinate frame and Fourier transformed in time (ω) and azimuth (m). The resulting equations are discretized as discussed in § 2.

The eddy-viscosity model we use, discussed in § 2, simply replaces µ in equations (B 4) to (B 6) with µ + µ T (x, r).

Appendix C. Optimizing in an input and output framework

In resolvent analysis, it is often useful to restrict the input and output spaces by writing

where C transforms the state vector to a desired output space y and B maps a smaller dimensional input space to the state space. Here we show that such additions do not hinder the generality of the optimization presented in this manuscript.

The structure of the cost function does not change,

but the SPOD modes, ψ, and the resolvent modes, u, are now computed considering the observable y, and the appropriate norms for the input and output space are defined by including weighting matrices W y andW f , respectively. The Lagrangian functional also takes a similar form as § 3.1,

, where, ũ1 , ṽ1 , σ1 are the Lagrange multipliers. The effective composition of the functional is identical to that of the full-state optimization as it is composed of the cost function, the forward solution, the resolvent eigenvalue problem, and a normalization constraint. Taking variations with respect to each variable, with exception to the eddy-viscosity term, results in the following system of equations,