# Eduction of coherent structures from schlieren images of twin jets using SPOD informed with momentum potential theory in the spectral domain

Iván Padilla-Montero, Daniel Rodríguez, Vincent Jaunet, Peter Jordan

► To cite this version:

Iván Padilla-Montero, Daniel Rodríguez, Vincent Jaunet, Peter Jordan. Eduction of coherent structures from schlieren images of twin jets using SPOD informed with momentum potential theory in the spectral domain. Theoretical and Computational Fluid Dynamics, 2024, Special Issue: Data-driven and Physics-based Modelling of Coherent Structures in Turbulent Shear Flow, 38 (3), pp.375-401. 10.1007/s00162-024-00699-w. hal-05374772

HAL Id: hal-05374772

https://hal.science/hal-05374772v1

Submitted on 20 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Iván Padilla-Montero · Daniel Rodríguez · Vincent

Jaunet · Peter Jordan

Eduction of coherent structures from schlieren
images of twin jets using SPOD informed with
momentum potential theory in the spectral domain

Received: date / Accepted: date

**Abstract** This work presents a methodology to extract coherent structures from high-speed schlieren images of turbulent twin jets which are more physically interpretable than those obtained with currently existing techniques. Recently, Prasad & Gaitonde [44] introduced an approach which employs the momentum potential theory of Doak [15,28] to compute potential (acoustic and thermal) energy fluctuations from the schlieren images by solving a Poisson equation, and combines it with spectral proper orthogonal decomposition (SPOD) to educe coherent structures from the momentum potential field instead of the original schlieren field. While the latter field is dominated by a broad range of vortical fluctuations in the turbulent mixing region of unheated high-speed jets, the momentum potential field is governed by fluctuations which are intimately related to acoustic emission, and its spatial structure in the frequency domain is very organized. The proposed methodology in this paper improves the technique of [44] in three new ways. First, the solution of the Poisson equation is carried out in the frequency-wavenumber domain instead of the time-space domain, which simplifies and integrates the solution of the Poisson equation within the SPOD framework based on momentum potential fluctuations. Second, the issue of solving the Poisson equation on a finite domain with *ad hoc* boundary conditions is explicitly addressed, identifying and removing those unphysical harmonic components introduced in the solution process. Third, the solution of the SPOD problem in terms of momentum potential fluctuations is used to reconstruct schlieren SPOD fields associated with each mode, allowing the visualization of the obtained coherent structures also in terms of the density gradient. The method is applied here to schlieren images of a twin-jet configuration with a small jet separation at two supersonic operation conditions: a perfectly-expanded and an overexpanded one. The SPOD modes based on momentum potential fluctuations retain the wavepacket structure including the direct Mach-wave radiation, together with upstream- and downstream-traveling acoustic waves, similar to SPOD modes based on the schlieren images. However, for the same dataset, they result in a lower-rank decomposition than schlieren-based SPOD and provide an effective separation of twin-jet fluctuations into independent toroidal and flapping oscillations that are recovered as different SPOD modes. These coherent structures are more consistent with twin-jet wavepacket models available in the literature than those originally obtained with direct schlieren-based SPOD, facilitating their interpretation and comparison against theoretical analyses.

**Keywords** Twin jets · SPOD · Acoustics · Schlieren · Momentum potential theory

Iván Padilla-Montero · Daniel Rodríguez

ETSIAE-UPM, School of Aeronautics, Universidad Politécnica de Madrid, Pza. Cardenal Cisneros 3, 28040

Madrid, Spain

E-mail: ivan.padilla@upm.es

Vincent Jaunet · Peter Jordan

Département Fluides, Thermique et Combustion, Institut Pprime, CNRS - Université de Poitiers - ISAE-ENSMA, 86036 Poitiers, France
# 1 Introduction

Noise pollution remains one of the major environmental health concerns associated with aviation, owing to its significant impact on the activity of human and animal life. Jet-engine noise is one of the main contributors to this problem, affecting both urban environments with civil aviation operations as well as military environments such as launch centers or aircraft carriers. Supersonic multi-jet engines are frequently encountered in the propulsion systems of rockets and modern high-speed aircraft. Among these, twin-jet configurations are the most common. Closely-spaced jets are known to strongly interact at the hydrodynamic and acoustic levels, giving rise to more complex flow structures than single round jets and distinct noise radiation patterns [4,30].

The link between the far-field sound radiated by high-speed jets and the turbulent fluctuations found in their core and mixing regions was soon recognized as a crucial point for the understanding of sound-generation mechanisms, and has been the subject of considerable research [27]. Multiple studies on single round jets have shown that the radiated sound is highly directional for both subsonic and perfectly-expanded supersonic jets [13,65,6], and that it is associated with large-scale, low-frequency fluctuations found in the mixing region [29,8]. Such large-scale fluctuations, today known as coherent structures, were first observed by Crow & Champagne [14]. These structures resemble the instability waves that develop in harmonically-forced jets, prompting their modeling by means of linear stability theory [12,36] and introducing the term “wavepackets” to refer to them.

During the last two decades, a successful modeling of wavepackets in isolated round jets was achieved by means of linear stability calculations, and their relevance in the coherent dynamics was further demonstrated through experimental comparisons and high-fidelity simulations [64,25,7]. Vortex sheet/finite thickness models [66], parabolized stability equations (PSE) [43, 47, 50, 58] and, more recently, one-way Navier-Stokes equations [68] and resolvent analysis [23, 26, 54] have been established as valid wavepacket modeling strategies. These studies have led to big advancements in the understanding of sound generation mechanisms through the modeling of Kelvin-Helmholtz instabilities and of acoustic resonances involving duct modes and shock-cell interactions in supersonic jets [67,18].

Analogous wavepacket models for subsonic and supersonic twin jets have only been developed in the last few years, as they pose additional challenges in terms of complexity and computational cost compared to single round jets: the azimuthal Fourier decomposition of the flow field is no longer possible due to the loss of axial symmetry of the flow structure of each jet, thus requiring the use of three-dimensional techniques. Studies based on local cross-plane linear stability theory with two inhomogeneous directions [49, 40, 51], along with plane-marching PSE [48,42] and twin-jet vortex sheet/finite thickness models [57, 38, 16, 63] have been used to characterize the instabilities that govern the twin-jet system and model the corresponding wavepackets.

Despite the achieved progress, current wavepacket models are not yet able to provide a complete description of the mechanics of sound generation in twin-jet systems. Linear models present limitations when mechanisms involving multiple interacting waves or coupling between the fluctuations of each jet are present. The interpretation of results is further complicated in configurations with small jet separation as the wavepackets become heavily deformed, departing from axial symmetry and complicating the identification of the oscillation modes. Therefore, hand-in-hand experimental investigations are necessary to validate the models as well as to characterize and provide physical understanding on mechanisms that have not yet been described computationally.

With the increasing ability to record and process high-resolution, high-speed experimental measurements, data-driven techniques introduced some decades ago have recently evolved and gained attention due to their potential to extract information from experimental data. In particular, the spectral proper orthogonal decomposition (SPOD) [3,11,69] offers a powerful framework to extract coherent structures from time-resolved experimental visualizations, such as those obtained from high-speed schlieren imaging. Recent investigations have employed this technique to study the screech resonance mechanism in single round jets with very satisfactory results [18,34]. For twin-jet configurations, the application of SPOD to experimental schlieren datasets has also been successful to describe screech resonances [20, 39, 45, 73, 41], where the associated coherent structures are composed of waves which, within each jet, feature a spatial structure analogous to a single azimuthal mode in the equivalent isolated jet. Conversely, at other frequencies where mixing-layer noise is the dominant mechanism, SPOD based on schlieren images is not well suited for separating the toroidal and flapping modes of oscillation of the
<figure><img src="image_2.png" /><figcaption>Fig. 1 Twin jet geometry and associated parameters.</figcaption></figure>

twin-jet system. As will be shown in this study, both types of oscillation are not clearly distinguishable in the structure of the leading SPOD modes corresponding to mixing-noise frequencies.

In a recent work, Prasad & Gaitonde [44] presented an approach to extract coherent structures from high-speed schlieren images using a quantity derived from the momentum potential theory formulation of Doak [15]. Jordan et al. [28] were the first to apply the momentum potential formulation to a complex model flow problem, consisting of a wavepacket with solenoidal and irrotational components, to investigate the source and flux terms involved in Doak's fluctuation-energy balance and their role in the generation of sound radiation. Later on, Unnikrishnan et al. [71, 72] applied the technique to a supersonic cold jet flow field obtained by means of large-eddy simulation, while Prasad & Morris [46] used it to study the performance of fluid inserts for the control of heated supersonic jets. In the approach of [44], the Helmholtz decomposition of the momentum density is used to derive a Poisson equation relating the schlieren fluctuation field with the streamwise gradient of the momentum potential fluctuation field integrated along the line of sight. By applying SPOD to the derived momentum potential fluctuations, coherent structures representing the acoustic and thermal energy components of the flow may be obtained. In [44], this technique is applied to schlieren visualizations of single round jets and twin rectangular jets, inducing coherent structures associated with sound generation and with the waves responsible for resonant feedback loops. Nevertheless, the performance and limitations of this approach compared to the direct use of schlieren images to extract coherent information have not yet been characterized. In particular, the impact of imposing *ad hoc* boundary conditions for the solution of the Poisson equation in a truncated domain where the flow field at the boundaries is not purely irrotational remains unclear [55, 56].

In this work, the idea introduced by Prasad & Gaitonde [44] is adopted to extract coherent structures from high-speed schlieren measurements in twin supersonic jets with an enhanced interpretability. A new methodology is proposed to advance the technique of [44] in three aspects: (i) integrating the calculation of the streamwise gradient of the momentum potential field within the SPOD algorithm by solving the Poisson equation in the spectral (frequency-wavenumber) domain; (ii) paying special attention to the influence of the *ad hoc* boundary conditions imposed on the Poisson equation, providing a filtering strategy to eliminate unphysical harmonic waves introduced in the process; (iii) reconstructing schlieren SPOD modes using the eigeninformation of the momentum potential SPOD problem. The newly obtained coherent structures are then compared against the SPOD structures provided by the direct use of schlieren images, allowing a quantification of the benefits attained from the use of the derived momentum potential fluctuations for this problem.

The remainder of the paper is organized as follows. Section 2 presents the experimental twin-jet setup employed to obtain schlieren visualizations. Section 3 describes the methodology used for the calculation of coherent structures, outlining the derivation of the Poisson equation relating schlieren and momentum potential fluctuations, the SPOD algorithm and the integration of the calculation of the momentum potential field within it. Section 4 presents SPOD results for two different twin-jet operating conditions and highlights the differences between the coherent information obtained from momentum potential fluctuations against schlieren fluctuations. Finally, concluding remarks are provided in § 5.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Experimental setup: (left) overall view of the twin-jet system; (middle) close-up view of the twin nozzles; (right) some elements of the schlieren system: flat mirror, knife edge, lens and camera.</figcaption></figure>

## 2 Experimental setup

The studied twin-jet configuration is represented in figure 1. The jets are generated by two identical round convergent-divergent nozzles. The nozzle geometry has been designed at Institut Pprime (CNRS-Université de Poitiers-ISAE-ENSMA), and follows a truncated ideal contour (TIC) profile with an exit diameter of $D = 0.025$ m and an exit-to-throat area ratio of $A_e/A_t = 1.225$. The center of each of the nozzles is located along the y-axis ($z = 0$) and the exit of each of them is located at $x = 0$. The spacing between the axis of each nozzle is denoted by $s$.

The experimental schlieren visualizations have been performed at the PROMÉTÉE platform of Institut Pprime (CNRS-Université de Poitiers-ISAE-ENSMA). The facility employed is the T200 compressible wind tunnel, which is powered by a 200 bar compressed air network and can reach operational conditions up to an isentropic Mach number $M_j = 2$ for the employed nozzle geometry. A heating system based on a series of tanks with heated nickel balls is used to increase and maintain the total temperature of the air arriving to the nozzles. The first and second pictures shown in figure 2 display the twin-jet experimental system built in the facility, which is placed inside a semi-anechoic room. The jet spacing considered for this study is $s/D = 1.76$, which is the minimum possible nozzle spacing according to the outer surface dimensions of the designed nozzles.

The jets are subject to different nozzle pressure ratios (NPR) with a fixed total temperature of $T_{t0} = 300$ K. The pressure ratio considered, i.e. the ratio of the total pressure in the reservoir $p_{t0}$ to the ambient pressure $p_\infty$ is defined in terms of the isentropic jet Mach number $M_j$ through the isentropic relation $p_{t0}/p_\infty = (1 + 0.5(\gamma - 1)M_j^2)^{(\gamma-1)/\gamma}$. Since the jets are not heated, the flow acceleration within the nozzle results in jet static temperatures lower than the ambient temperature. The operation condition closest to the perfectly-expanded regime that could be achieved in practice is obtained at $M_j = 1.54$, calibrated experimentally by means of flow visualization. In this work, two different operating conditions are considered, namely $M_j = 1.54$ and $M_j = 1.26$, the latter corresponding to an overexpanded regime. The corresponding Reynolds number based on the nozzle exit diameter and the jet exit flow conditions for $M_j = 1.54$ is $Re = 1.40 \times 10^6$.

The optical system consists of a classical Z-type schlieren setup. Some of its elements are shown in the third picture of figure 2. A continuous light source is provided by a 60 W LED, which goes through an aperture that prevents any light from the source to directly enter the test section. Two parabolic mirrors with a diameter of 30 cm and a focal length of 3 m are used to produce a collimated light beam which, through two circular apertures, traverses the test section along the z direction according to the reference frame shown in figure 1. Two additional flat mirrors with a diameter of 12 cm each are used to accommodate the optical path of the system into a Z shape, allowing for a more compact experimental setup. A vertical knife edge is placed at the focal length of the second parabolic mirror; in consequence, the light intensity field in the recovered schlieren images is proportional to the streamwise density gradients in the flow $(\partial\rho/\partial x)$.

The camera employed to record the images at high-speed is a Phantom v2640. For each operating condition, 30,000 snapshots are recorded at a sampling frequency of $f_s = 68$ kHz with a spatial
<figure><img src="image_4.png" /><figcaption>Fig. 3 Schlieren visualizations of the twin jet at $M_j = 1.26$ (overexpanded regime): (left) sample instantaneous snapshot; (right) mean field obtained by averaging 30,000 snapshots.</figcaption></figure>

resolution of 352 × 512 pixels. An exposure time 0.8 µs is used, which is small enough to freeze the convected flow disturbances while ensuring enough contrast in the image. The recorded images consist of an array of pixel intensity values ranging between 0 and 4095 (12-bit color depth).

The image's spatial window approximately comprises an area of 3D away from each nozzle axis along the y direction and 10D downstream of the nozzle exit along the x direction. This size was found to offer a good compromise between spatial resolution and window size for the current setup, resulting in a pixel size of 0.02D. Figure 3 displays an instantaneous schlieren snapshot obtained at $M_j = 1.26$, as well as the corresponding mean field obtained by averaging all 30,000 snapshots.

## 3 Methodology

In the following, instantaneous flow quantities $q$ are split into a time-stationary mean component and a time-dependent fluctuating component: $q(x, y, z, t) = \bar{q}(x, y, z) + q'(x, y, z, t)$.

### 3.1 Doak's momentum potential theory

The momentum potential formulation proposed by Doak [15] is based on a Helmholtz decomposition of the form:

$$ \rho u_i = B_i - \frac{\partial \psi}{\partial x_i}, \qquad (1) $$

where $\rho u_i$ denotes the momentum density of the flow (density $\rho$ times velocity $u_i$), $B_i$ is the solenoidal momentum component ($\partial B_i / \partial x_i = 0$) and $\psi$ denotes the scalar potential associated with the irrotational momentum component ($\nabla \times \partial \psi / \partial x_i = 0$).

Introducing the decomposition (1) into the continuity equation

$$ \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u_i)}{\partial x_i} = 0, \qquad (2) $$

assuming a statistically time-stationary flow ($\partial \bar{\rho} / \partial t = 0$), and noting that $\rho(x_i, t) = \bar{\rho}(x_i) + \rho'(x_i, t)$ and $\psi(x_i, t) = \bar{\psi}(x_i) + \psi'(x_i, t)$, yields, on one hand:

$$ \frac{\partial^2 \psi'}{\partial x_i^2} = \frac{\partial \rho'}{\partial t}, \qquad (3) $$
which is a Poisson equation relating the two scalar fluctuation fields $\rho'$ and $\psi'$, and, on the other hand:

$$ \frac{\partial^2 \bar{\psi}}{\partial x_i^2} = 0, \qquad (4) $$

which implies that the scalar potential field has zero mean, i.e. $\psi(x_i, t) = \psi'(x_i, t)$ (see [15]).

Helmholtz's decomposition theorem expressed in the form of equation (1) holds for vector fields that decay to zero sufficiently fast at infinity [2,62]. This requirement is usually not satisfied in the finite domains employed for the numerical simulation or experimental visualization of flow fields of practical interest, such as the twin-jet configuration under study. Certainly, the momentum density of the twin-jet flow field does not vanish at the boundaries of the domain shown in figure 3. In these cases, the truncation of the domain introduces a third component in the decomposition which corresponds to an harmonic solution component (see Theorem 2 in [56]):

$$ \rho u_i = B_i - \frac{\partial \psi}{\partial x_i} - H_i. \qquad (5) $$

The harmonic component, denoted by $-H_i$, is both solenoidal and irrotational [2]. As a result, it also has an associated scalar potential $H_i = \partial\psi_h/\partial x_i$ that satisfies Laplace's equation $(\nabla \cdot H_i = \partial^2\psi_h/\partial x_i^2 = 0)$ in the domain of analysis. This component represents the impact of the truncation of the domain on the Helmholtz decomposition [55,56]. Consequently, as the schlieren images are contained in a truncated domain with boundaries at which the flow momentum does not vanish, $H_i$ will also be recovered as part of the solution of equation (3). In those cases where ad hoc boundary conditions are enforced in the solution of the Poisson equation, it is crucial to pay attention to the impact that the harmonic field has on the physically-relevant part of the solution, as discussed in section 3.4.

## 3.2 Relation between schlieren fluctuations and scalar momentum potential fluctuations

According to the experimental set-up described in section 2, the recorded schlieren images represent streamwise density gradients integrated along the line of sight (z axis), that is

$$ \sigma(x, y) = \int \frac{\partial \rho(x, y, z)}{\partial x} dz. \qquad (6) $$

Following Prasad & Gaitonde [44], equation (3), owing to its linear nature, can be differentiated with respect to x to obtain:

$$ \frac{\partial^2}{\partial x_i^2} \left( \frac{\partial \psi'}{\partial x} \right) = \frac{\partial}{\partial t} \left( \frac{\partial \rho'}{\partial x} \right). \qquad (7) $$

If integrated along the schlieren line of sight, equation (7) directly relates the fluctuation of the schlieren field with the streamwise derivative of the potential fluctuation field $\psi'$ integrated along the line of sight, that is

$$ \frac{\partial^2 \Theta'}{\partial x_i^2} = \frac{\partial \sigma'}{\partial t}, \qquad (8) $$

where $\sigma' = \int (\partial\rho'/\partial x) dz$ denotes the schlieren fluctuation field and $\Theta' = \int (\partial\psi'/\partial x) dz$.

As shown by Doak [15], the fluctuating potential $\psi'$ can be written as the sum of acoustic and thermal components, $\psi' = \psi'_A + \psi'_T$, which are directly related with the pressure and the entropy fluctuation fields, respectively. This decomposition bears strong resemblance with the generalized potential disturbance energy definition introduced by Chu [10], which consists of an acoustic component associated with compression due to pressure fluctuations, and a thermal component associated with heat exchange due to entropy spottiness. In this regard, $\psi'$, and by extension $\Theta'$, can be interpreted as
quantities describing the potential fluctuation energy in the system. Previous works [71, 72, 46] indicate that the thermal component is non-radiating in cold supersonic turbulent jets. As a result, outside of the jet core, $\Theta'$ represents the acoustic energy present in such flows.

Equation (8) may be interpreted as a means to calculate $\Theta'$ from experimental schlieren measurements. However, for this statement to be rigorous, the boundary conditions imposed on equation (8) must satisfy the irrotationality condition inherent in $\partial\psi'/\partial x_i$. In other words, the calculation of the boundary values of $\Theta'$ requires the integration of equation (1) along far-field boundaries at which the momentum density of the flow is purely irrotational ($\rho u_i = -\partial\psi'/\partial x_i$). In general, this is not possible when dealing with experimental datasets, as measurements of the momentum density are not available. For numerical datasets, although the momentum density field is known, the domain boundaries rarely contain a purely irrotational flow. In convective flow simulations, and particularly in jet calculations, the shear layer usually exits the domain through the outflow boundary, where the flow field is by no means irrotational. Therefore, the aforementioned condition can only be rigorously satisfied in certain particular cases such as the solenoidal wavepacket investigated by Jordan et al. [28]. Despite this difficulty, the solution of the Poisson equation can still be carried out in practice using an ad hoc treatment at the domain boundaries, at the expense of introducing unphysical waves in the solution through the harmonic component (see section 3.1). This is, to the best of the authors' knowledge, the approach adopted by [71, 72, 46, 44, 45].

In this work, the proposed methodology enables the identification of the harmonic solution component that arises from the domain truncation and associated boundary conditions, and allows to remove its influence from the physically-relevant solution component, as described in section 3.4.

### 3.3 Spectral proper orthogonal decomposition

The spectral proper orthogonal decomposition (SPOD) algorithm as described by Towne et al. [69] is employed here to extract the spatio-temporal coherent structures from the time-resolved experimental data. The algorithm employs Welch's method to average the spectral information over multiple realizations of the flow and as a result obtain convergent estimates of the cross-spectral density tensor. In this study, the number of realizations in use is much smaller than the number of grid points (image pixels) of the schlieren datasets, and therefore, for a given frequency, the snapshot method [59] variant of the SPOD eigenvalue problem is used:

$$ \tilde{\mathbf{Q}}_k^H \mathbf{W} \tilde{\mathbf{Q}}_k \mathbf{\Phi}_k = \mathbf{\Phi}_k \mathbf{\Lambda}_k, \qquad (9) $$

where $\tilde{\mathbf{Q}}_k$ is the matrix whose columns contain the temporal discrete Fourier transform (DFT) of each segment of the total time-series used in Welch's averaging (considered as an independent realization) at the $k$th frequency, $\mathbf{W}$ is the matrix that contains the weight and the numerical quadrature defining the SPOD spatial inner product, $\mathbf{\Phi}_k$ is a matrix whose columns contain the eigenvectors of the problem for the $k$th frequency and $\mathbf{\Lambda}_k = \text{diag}(\lambda_1, \lambda_2, \dots)$ is a diagonal matrix containing the associated eigenvalues for that particular frequency. The superscript $H$ denotes the complex-conjugate transpose. The eigenvalues found in $\mathbf{\Lambda}_k$ directly represent the power spectral density associated with each of the SPOD modes. The SPOD modes for the $k$th frequency, denoted by $\mathbf{\Psi}_k$, are obtained by means of the following relation:

$$ \mathbf{\Psi}_k = \tilde{\mathbf{Q}}_k \mathbf{\Phi}_k \mathbf{\Lambda}_k^{-1/2}. \qquad (10) $$

The assembly of the matrix $\tilde{\mathbf{Q}}_k$ for each frequency requires prior calculation of the temporal DFT of each realization of the flow. In this case, each realization consists of a data matrix (block) with a segment of the time-series of schlieren fluctuation ($\sigma'$) or $\Theta'$ snapshots of equal length. In this work, each dataset consists of 30,000 instantaneous snapshots, which are divided into 57 blocks of 1024 snapshots each with a 50% overlap. Each block is windowed using a Hamming window to reduce spectral leakage.

Given the nature of the experimental datasets (pixel intensities), the definition of an integral energy norm is not meaningful. Instead, a weighted 2-norm based on a trapezoidal integration rule (see [53]) is employed to serve as a numerical quadrature for the SPOD inner product.
This work focuses on those coherent structures that describe the spatio-temporal dynamics of the twin jet system when both jets are coupled, that is, structures that involve both jets. For this purpose, the symmetry of the system is exploited to generate symmetric and antisymmetric fluctuation fields. When schlieren realizations are considered, each snapshot is split by the line at $y = 0$ to create two new datasets, namely, a symmetric dataset $\sigma_s = (\sigma_u + \sigma_l)/2$ and an antisymmetric one $\sigma_a = (\sigma_u - \sigma_l)/2$, with $\sigma_u$ and $\sigma_l$ respectively denoting the upper ($y > 0$) and lower ($y < 0$) halves of the original schlieren images. The SPOD is then applied to both the symmetric and antisymmetric datasets separately.

In the results presented in this work, the SPOD problem (9) is solved both in terms of the cross-spectral density of $\sigma'$ or in terms of the cross-spectral density of $\Theta'$. This is done by constructing the $\tilde{Q}_k$ matrix using either the temporal DFT of schlieren realizations or the temporal DFT of $\Theta'$ realizations. The DFT of schlieren realizations is straightforward to obtain by means of the fast Fourier transform (FFT). The approach employed to obtain the DFT of $\Theta'$ fields for each realization is described in the following section.

### 3.4 Calculation of $\Theta'$ and its integration within the SPOD framework

Poisson equation (8) can be solved directly in the space-time domain by employing discretizations of the temporal and spatial derivatives and solving a linear system of equations, as done by e.g. Prasad & Gaitonde [44]. In this work, an alternative procedure is employed which solves the equation in the frequency-wavenumber domain by performing discrete Fourier transforms (DFT) in time and on the Cartesian two-dimensional domain of the schlieren images. The reasons for this choice are elaborated next.

First, solving for $\Theta'$ in the wavenumber domain prevents discretizing the spatial derivatives of the Laplacian operator and the consequent solution of a linear system of equations, in exchange of performing the FFT along $x$ and $y$. Second, working in the frequency domain simplifies the calculation of the time derivative of the schlieren fluctuations with high accuracy [37]. In the time domain, accurate time-derivative calculations require the use of very high-order differentiation schemes. As shown in Appendix B, computing the time derivative of $\sigma'$ by means of a low-order finite difference scheme results in a significant damping of the high-frequency content present in the schlieren snapshots, which slowly decreases as the order of the scheme increases. Standard central finite difference schemes, however, have accuracy limitations when the order becomes very high, as the finite difference stencils become ill-conditioned [52]. Therefore, more sophisticated high-order differentiation methods (e.g. Padé finite difference schemes [37] or collocation methods [70]) are required for accurate time-derivative calculations in the temporal domain, which carry an associated increase in computational complexity. Third, since the use of $\Theta'$ realizations in the SPOD algorithm requires the temporal Fourier transform of the $\Theta'$ field to construct the $\tilde{Q}_k$ matrix (see section 3.3), obtaining $\Theta'$ directly in the frequency domain does not involve additional cost, as it only brings forward the necessary temporal DFT. Appendix C provides a comparison of the computational cost of the SPOD solution when calculating $\Theta'$ realizations in the space-time domain, against the cost of the SPOD algorithm computing $\Theta'$ realizations in the frequency-wavenumber domain following the procedure detailed below.

Performing the temporal DFT, equation (8) becomes:

$$ \frac{\partial^2 \tilde{\Theta}'}{\partial x_i^2} = i\omega \tilde{\sigma}', \qquad (11) $$

where $\omega = 2\pi f$ is the angular frequency and the tilde symbol denotes Fourier-transformed quantities in time, e.g. $\tilde{q}'(x, y, \omega)$. Performing the DFT in both $x$ and $y$, the resulting form of the equation in the frequency-wavenumber domain is

$$ -(k_x^2 + k_y^2) \hat{\Theta}' = i\omega \hat{\sigma}', \quad \text{or} \quad \hat{\Theta}' = -\frac{i\omega \hat{\sigma}'}{k_x^2 + k_y^2}, \qquad (12) $$

where $k_x$ and $k_y$ respectively denote the wavenumbers along $x$ and $y$ and the hat symbol refers to Fourier-transformed quantities in time and space, $\hat{q}'(k_x, k_y, \omega)$. Equation (12) provides an algebraic
expression to directly compute $\hat{\Theta}'$ for each frequency and wavenumber pair. In order to obtain $\hat{\Theta}'$, which is the quantity of interest for the SPOD algorithm, the inverse DFT of $\hat{\Theta}'$ in space is carried out.

The aforementioned procedure imposes periodicity of the solutions at the boundaries of the schlieren domain. In order to artificially place the periodic boundaries further away from the actual boundaries dictated by the schlieren image size, the temporal Fourier-transformed schlieren fluctuation $\tilde{\sigma}'$ is zero-padded along the $x$ and $y$ directions before transforming it to the wavenumber domain. To prevent a discontinuous jump in $\tilde{\sigma}'$ when extending its spatial domain with zeros, a spatial window can be applied to the $\tilde{\sigma}'$ field prior to zero padding, thus enforcing a smooth transition of the field towards zero at the domain boundaries. Here, a Planck-taper window [35] is considered for that purpose, which enables a rapid but smooth transition to zero at the domain boundaries without significantly altering the interior values. The values of the window function $w$ for a one-dimensional signal of $N + 1$ points can be expressed as follows:

$$ \left\{ \begin{array}{ll} w_0 = 0, & \\ w_j = \left[ 1 + \exp\left(\frac{\varepsilon N}{j} - \frac{\varepsilon N}{\varepsilon N - j}\right) \right]^{-1}, & \text{if } 1 \le j < \varepsilon N \\ w_j = 1, & \text{if } \varepsilon N \le j \le N/2 \\ w_j = w_{N-j}, & \text{if } N/2 < j \le N \end{array} \right. \qquad (13) $$

where $0 < \varepsilon \le 0.5$ is the so-called tapering parameter. To apply this bump window to the two-dimensional $\tilde{\sigma}'$ fields, a two-dimensional Planck-taper window must be generated. This is achieved by first generating two one-dimensional Planck-taper windows, one with $N_x$ points and another with $N_y$ points, with $N_x$ and $N_y$ being the number of points of the original schlieren images along $x$ and $y$, respectively. Then, each of these functions are replicated $N_y$ and $N_x$ times, respectively, to generate two two-dimensional functions which are multiplied together to obtain a two-dimensional bump function that can be applied to each $\tilde{\sigma}'$ field directly. The sensitivity of the SPOD solution to the tapering parameter $\varepsilon$ is small for the analyzed datasets. It is found to be stronger for low frequencies, which are associated with structures of larger spatial wavelength.

After computing $\hat{\Theta}'$ in the frequency-wavenumber domain and inverting the spatial DFT along $x$ and $y$ to obtain $\hat{\Theta}'$, the spatial domain is reduced back to the original one by removing the spatial points corresponding to the previously added zero padding.

As introduced before, the calculation of $\hat{\Theta}'$ from schlieren snapshots by means of the described methodology introduces harmonic unphysical solution components in the obtained $\hat{\Theta}'$ fields due to the violation of the irrotationality condition at the boundaries of the schlieren domain. The harmonic unphysical component is associated with the imposition of ad hoc boundary conditions (in this case, periodic boundary conditions) in the artificially truncated physical domain; it is found to be sensitive to changes in the domain size and its energy to be concentrated in small streamwise wavenumbers, featuring streamwise phase speeds $c_{ph} = \omega/k_x$ that are highly supersonic and do not represent physically-sound components of the twin-jet system. Appendix A illustrates the signatures of these unphysical harmonic waves, and how they can be clearly distinguished from physically-meaningful components. No unphysical energetic components have been observed in the $y$ direction, which is attributed to the fact that for this problem, the main rotational flow regions at the boundaries are found at the streamwise boundaries (mainly the downstream boundary of the schlieren images).

In order to prevent these harmonic components from contaminating the SPOD modes to be computed subsequently, the energy of those $\hat{\Theta}'$ components with a streamwise wavenumber $k_x$ that is smaller than a given cut-off value is set to zero. Such cut-off value is frequency dependent, and is fixed by specifying a supersonic streamwise phase speed $c_{ph,c}$ threshold. This filtering procedure is applied in the wavenumber domain, just before inverting the spatial DFT to obtain a $\hat{\Theta}'$ fluctuation field that does not contain the energetic part of the unphysical harmonic component.

The proposed approach to construct the SPOD eigenvalue problem (9) based on $\Theta'$ can be summarized in the following steps:

1. Compute the temporal FFT of schlieren fluctuation snapshots ($\tilde{\sigma}'$) for each realization (block).

2. For each frequency and each block:
(a) Add a zero-padding along $x$ and $y$ to the $\tilde{\sigma}'$ field. To prevent a discontinuous jump in $\tilde{\sigma}'$ when padding, a spatial window should be applied beforehand. In this work, a Planck-taper window with $\epsilon = 0.1$ has been employed.

(b) Compute the spatial FFT of the zero-padded $\tilde{\sigma}'$ field along $x$ and $y$.

(c) Calculate the $\hat{\Theta}'$ field using equation (12).

(d) Remove the most energetic unphysical harmonic waves by setting to zero the region of the $\hat{\Theta}'$ field comprised between $-\omega/c_{ph,c} < k_x < \omega/c_{ph,c}$, where $c_{ph,c}$ is the chosen phase-speed threshold. In this study, a value of $c_{ph,c} = 1.2c_\infty$ has been used, where $c_\infty$ denotes the freestream speed of sound. The sensitivity of the SPOD results to the choice of $c_{ph,c}$ is reported in Appendix A.

(e) Compute the inverse spatial FFT of the filtered $\hat{\Theta}'$ field to obtain $\tilde{\Theta}'$.

(f) Reduce the spatial window of the $\hat{\Theta}'$ field back to the original size of the schlieren snapshots.

(g) If desired, create symmetric and antisymmetric $\hat{\Theta}'$ fields as $\hat{\Theta}'_s = (\hat{\Theta}'_u + \hat{\Theta}'_l)/2$ and $\hat{\Theta}'_a = (\hat{\Theta}'_u - \hat{\Theta}'_l)/2$.

3. For each frequency, build the $\tilde{Q}_k$ matrix using the computed $\tilde{\Theta}'$ fields for all realizations.

### 3.5 Calculation of schlieren SPOD modes for the $\Theta'$ cross-spectral density

The snapshot formulation of POD allows the computation of flow-field variables, corresponding to the SPOD modes, that are not involved in the definition of the cross-spectral density matrix [21, 5, 58, 61, 60, 31, 33, 32]. In particular, it allows the computation of the schlieren fluctuation field $\sigma'$ (the line-of-sight integrated $\partial\rho'/\partial x$) associated with the SPOD modes corresponding to the cross-spectral density of $\Theta'$.

This is accomplished using the eigenvectors and eigenvalues of the $\Theta'$ decomposition ($\Phi_{k,\Theta'}$ and $\Lambda_{k,\Theta'}$, respectively) combined with the matrix of frequency-domain realizations for the schlieren fields ($\tilde{Q}_{k,\sigma'}$):

$$ \Psi_{k,\sigma'} = \tilde{Q}_{k,\sigma'} \Phi_{k,\Theta'} \Lambda_{k,\Theta'}^{-1/2}. \quad (14) $$

Note that this procedure is equivalent to solving an augmented SPOD eigenvalue problem consisting of a $\tilde{Q}_k$ matrix containing both schlieren and $\Theta'$ realizations, and then using a weight matrix that is non-zero only for the $\Theta'$ components.

## 4 Results

This section presents results obtained by applying the proposed methodology to the twin-jet schlieren visualizations acquired for two different operating conditions: $M_j = 1.54$, corresponding approximately to a perfectly-expanded jet condition, and $M_j = 1.26$, corresponding to an over-expanded jet condition. SPOD results for the two definitions of the cross-spectral density matrix are presented. SPOD results based directly on the schlieren images ($\sigma'$) are denoted as $C_{\sigma'}$, while those based on $\Theta'$ are denoted by $C_{\Theta'}$.

### 4.1 Perfectly-expanded condition

Perfectly-expanded jets do not present strong shock waves and their acoustic signature is characterized by broadband mixing-layer noise. In turn, schlieren images are expected to be dominated by a large range of vortical fluctuations at all frequencies, making their direct spectral analysis by means of DFT impractical. SPOD based on the schlieren images is nevertheless anticipated to recover the low-rank dynamics associated with those highly-energetic coherent $\sigma'$ fluctuations in the mixing-noise frequency band, which would manifest in the SPOD spectra as a bigger separation between the spectral density $\lambda_i$ of the leading SPOD modes and the spectral density of the rest of modes [54].
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Comparison of the SPOD spectra for $M_j = 1.54$. (left) $C_{\sigma'}$ SPOD; (right) $C_{\Theta'}$ SPOD; (top) symmetric case; (bottom) antisymmetric case. Each line corresponds to one SPOD mode. A grayscale is used to range from the most energetic mode to the least energetic one. The orange line represents the sum of the energy of all SPOD modes for each frequency.</figcaption></figure>

Figure 4 displays the symmetric and antisymmetric SPOD spectra obtained for $M_j = 1.54$ when using the correlations $C_{\sigma'}$ and $C_{\Theta'}$. The spectra are represented as a function of the Strouhal number, defined as $St = fD/u_j$, where $u_j$ is the isentropic jet velocity derived from $M_j$ and the isentropic nozzle exit temperature. For each frequency, the SPOD modes are ranked by their spectral density and are represented following a grayscale that ranges from the most relevant mode (mode 1, black) to the least relevant one (mode 57, white). In addition, a line representing the sum of the spectral density of all the SPOD modes ($\sum_i \lambda_i$) for each frequency is included in orange.

The SPOD spectra illustrate the expected broadband increase in spectral density ($St \approx 0.1$ to 1) associated with the coherent wavepackets that grow in the shear layers of the jets, and which, for the current case, dominate for an antisymmetric coupling of the jets. Both the schlieren and the $\Theta'$-based SPOD spectra show this trend, featuring the maximum spectral density at the same Strouhal numbers, $St = 0.3$ to 0.4. Nevertheless, two important differences between both spectra stand out:

(i) The sum of the spectral density for all modes (orange line) in the $C_{\Theta'}$ spectra is substantially reduced outside of the frequency range typical of mixing noise, in contrast to the $C_{\sigma'}$ spectra, where the total spectral density only decreases mildly outside of the mixing-noise range, featuring a broader band behaviour. This is consistent with the interpretation of $\Theta'$ as related to potential-energy fluctuations, and illustrates their dominance at the range of frequencies associated with mixing noise.

(ii) The $C_{\Theta'}$ SPOD results exhibit a lower-rank behavior than the schlieren-based counterpart, as indicated by the small difference between the orange line and the line corresponding to the leading SPOD mode (mode 1). Therefore, most of the spectral density for $C_{\Theta'}$ is accounted for by the leading SPOD mode.

The SPOD modes for both the $C_{\sigma'}$ and the $C_{\Theta'}$ SPOD solutions are illustrated in figures 5 and 6 for two different frequencies, respectively. Each figure shows the first (left column) and second (right column) SPOD modes in three different versions: (i) the schlieren field obtained from the schlieren SPOD ($C_{\sigma'}$) (first row); (ii) the $\Theta'$ field obtained from the $\Theta'$ SPOD problem ($C_{\Theta'}$) (second row); and (iii) the schlieren field reconstructed from the $\Theta'$ SPOD problem ($C_{\Theta'}$) as outlined in section 3.5 (third row). Due to the fact that the symmetric or antisymmetric fields are considered, only one of the two
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Contours of the real part of the antisymmetric SPOD modes for $M_j = 1.54$ at $St = 0.37$: (top) $\sigma'$ field, $C_{\sigma'}$ SPOD problem; (middle) $\Theta'$ field, $C_{\Theta'}$ SPOD; (bottom) $\sigma'$ field, $C_{\sigma'}$ SPOD. (left) mode 1; (right) mode 2. The grey regions represent the exterior nozzle surface and the part of the mirror support structure that is contained within the schlieren window. The black dashed lines depict the inner and outer nozzle lip lines for each jet.</figcaption></figure>

jets is shown in each contour plot. The quantity $\lambda_1/\lambda_2$ denotes the ratio between the energy of the first and the second SPOD modes, while $\lambda_1/\sum_i \lambda_i$ is the portion of the total energy contained in the first mode (and similarly for $\lambda_2/\lambda_3$ and $\lambda_2/\sum_i \lambda_i$).

Figure 5 shows the modes at the peak frequency of the broadband region for the antisymmetric coupling, located at $St = 0.37$. In agreement with previous investigations, the coherent structures resemble Kelvin-Helmholtz instability wavepackets that manifest as toroidal or flapping deformations of each jet column and present a noticeable downstream Mach wave radiation. Toroidal twin-jet modes are analogous to $m = 0$ modes in isolated round jets in the sense that they feature an axisymmetric amplitude function within each jet. Similarly, when a coupling exists between the oscillations of twin jets, it favors flapping motions over the $m = 1$ helical modes typical of single round jets [51]. The nature of the schlieren visualization prevents recovering oscillations that are anti-symmetric about the jet-containing plane; only symmetric (i.e. sinuous) or anti-symmetric (varicose) flapping motions can be observed. These oscillation modes are characterized by their vanishing amplitude along the axis of each jet.

However, owing to the mean flow interaction between both jets when they are closely spaced (such as in the configuration under study), the inner and outer shear layers of each jet have different thickness and shear magnitudes [24, 42]. As a result, the wavepackets no longer have a perfectly axisymmetric character, especially for lower frequencies which are associated with longer wavelengths, and the identification of the oscillation modes as $m = 0$ and $m = 1$ becomes less straightforward than for larger jet spacing [42]. In addition, the structures that characterize schlieren-based SPOD modes in twin jets may contain signatures of $m = 0$ and $m = 1$ waves mixed in the same mode [41]. This is the case of the schlieren SPOD modes shown in the top row of figure 5, where the supported coherent structures have an oblique orientation with respect to the jet axis and it is not evident whether their spatial structure corresponds to a toroidal or a flapping oscillation.
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Contours of the real part of the antisymmetric SPOD modes for $M_j = 1.54$ at $St = 0.61$: (top) $\sigma'$ field, $C_{\sigma'}$ SPOD problem; (middle) $\Theta'$ field, $C_{\Theta'}$ SPOD; (bottom) $\sigma'$ field, $C_{\sigma'}$ SPOD; (left) mode 1; (right) mode 2.</figcaption></figure>

Looking at the second row of figure 5, it can be observed that the structure of the $C_{\Theta'}$ SPOD modes also captures the Kelvin-Helmholtz waves and their associated acoustic radiation. Moreover, it reveals that the $C_{\Theta'}$ SPOD performs very well at separating the toroidal and the flapping structures into the first and second SPOD modes, respectively. This illustrates a key advantage of the methodology, namely that the $C_{\Theta'}$ SPOD can provide a more robust extraction of coherent information in the system than that based on $C_{\sigma'}$. Furthermore, note that $\lambda_1$ contains 74% of the total spectral density of the $C_{\Theta'}$ decomposition, while it only represents 20% of the total for $C_{\sigma'}$, which reflects the lower-rank behavior of the $C_{\Theta'}$ decomposition described before.

The effectiveness of the $C_{\Theta'}$ SPOD problem in separating and ranking the different coherent structures can then be brought to the schlieren field by reconstructing the schlieren modes corresponding to the $C_{\Theta'}$ SPOD modes, as presented in the third row of figure 5. The reconstructed schlieren field for the second SPOD mode shows a much clearer $m = 1$ structure than the original schlieren mode.

Figure 6 depicts the antisymmetric SPOD modes at a higher frequency ($St = 0.61$). In this case, it is more evident that $m = 0$ and $m = 1$ wavepackets are mixed in the original $C_{\sigma'}$ SPOD modes. The $C_{\Theta'}$ SPOD decomposition is found to effectively separate the structures into toroidal and flapping fluctuations in the first and second SPOD modes. The reconstructed schlieren SPOD fields represent an improved description of the coherent structures associated with sound generation in the twin-jet system.

## 4.2 Overexpanded condition

Next, an isentropic Mach number at which the jets are significantly overexpanded is considered ($M_j = 1.26$). For this condition, shock cells are present in the core of the jets (see figure 3), and the spectra is expected to be dominated by the screech resonance mechanism (see [17]) together with mixing noise.

Figure 7 shows the symmetric and antisymmetric SPOD spectra for this case, obtained solving both the $C_{\sigma'}$ and the $C_{\Theta'}$ SPOD problems. The highest energy in the spectra is now encountered in the
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Comparison of the SPOD spectra for $M_j = 1.26$. (left) $C_{\sigma'}$ SPOD; (right) $C_{\Theta'}$ SPOD; (top) symmetric case; (bottom) antisymmetric case.</figcaption></figure>

strong tonal peaks associated with the screech resonance. Three tonal peaks are visible, corresponding to the fundamental screech frequency ($St = 0.56$) and two of its harmonics. The tones are much more energetic for the symmetric spectra, indicating that the screech resonance for the condition under study favors a symmetric twin-jet coupling. Besides screech, the broadband signature associated with mixing noise can also be identified, once again showing higher spectral density for antisymmetric fluctuations. Except for the fundamental screech frequency, where both the $C_{\sigma'}$ and the $C_{\Theta'}$ spectra contain almost all of the energy in the first SPOD mode, the $C_{\Theta'}$ SPOD features a much lower-rank behavior also for this $M_j$ along the entire broadband region.

The SPOD modes for the over-expanded jet condition are presented in figures 8 to 12. First, the symmetric leading SPOD mode corresponding to the fundamental screech frequency ($St = 0.56$) is analyzed, illustrated in figure 8. As in the previous section, three different fields are represented: the schlieren field corresponding to $C_{\sigma'}$, the $\Theta'$ field corresponding to $C_{\Theta'}$, and the schlieren field corresponding to $C_{\Theta'}$.

The $C_{\sigma'}$ SPOD mode shows structures that correspond to toroidal wavepackets modulated by the shock cells, together with downstream and upstream acoustic radiation. Due to the strong resonance that originates the screech, in this case there is no difficulty in identifying the $m = 0$ structure in the amplitude function ($\lambda_1$ represents nearly 98% of the total energy). The $C_{\Theta'}$ SPOD mode structure also contains most of these features. Nonetheless, note that no signature of the shock-cell modulation is visible in the $\Theta'$ field, which illustrates how, within the frame of momentum potential theory, shock cells do not contribute directly to the irrotational fluctuation components.

In order to separate the upstream- and downstream-propagating wave components present in the SPOD mode associated to the fundamental screech frequency, the spatial DFT of the SPOD mode is calculated along $x$ to obtain its spatial structure as a function of the streamwise wavenumber $k_x$. This enables the separation of those energetic components with $k_x > 0$, i.e. positive phase speed, from those components with $k_x < 0$, i.e. negative phase speed. The spatial Fourier transform of the SPOD mode is then inverted by keeping only the $k_x > 0$ or the $k_x < 0$ components, resulting in the structures represented in figure 9. The right column shows waves that propagate downstream within the screech mode, consisting of a Kelvin-Helmholtz wavepacket and its associated acoustic radiation.
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Real part of the leading (mode 1) symmetric SPOD mode for $M_j = 1.26$ at $St = 0.56$, corresponding to the fundamental screech frequency: (top) $\sigma'$ field, $C_{\sigma'}$ SPOD problem; (middle) $\Theta'$ field, $C_{\Theta'}$ SPOD; (bottom) $\sigma'$ field, $C_{\Theta'}$ SPOD.</figcaption></figure>

The left column contains mostly waves that propagate upstream, and which take the form of acoustic waves with spatial support that extends into the freestream, the jet shear layers and in the subsonic regions of the jet core [19,34]. It is important to note that the actual direction of propagation of energy in the flow is dictated by the sign of the group velocity ($c_g = \partial\omega/\partial k_x$), rather than that of the phase velocity. Although most waves with negative phase speed also feature negative group velocity and thus propagate upstream (and vice-versa), there are also certain waves with negative phase speed that have a positive group velocity and move downstream, such as acoustic waves that travel within the supersonic jet core region.

The $C_{\Theta'}$ SPOD solution is successful in recovering both the negative and the positive phase-speed potential components of the screech resonance. In this case, the irrotational structures found in the $\Theta'$ mode follow a very similar spatial organization to the original schlieren SPOD mode, and as a result the reconstructed schlieren SPOD field is almost identical to the original one. Because its spatially organized structures are already well resolved in the original schlieren SPOD mode, the reconstructed field does not bring any advantage in terms of the interpretability of the coherent structures associated with the screech resonance.

The second frequency analyzed for this operating condition corresponds to the maximum spectral density in the mixing-noise broadband region, namely, $St = 0.36$, which is higher for the antisymmetric case. Figure 10 depicts the respective SPOD modes for this frequency. Similarly to the previous section, both the first (left column) and second (right column) SPOD modes are shown. In this case, the $C_{\sigma'}$ SPOD modes show apparently a combination of three structures with two different wavenumbers: one of them non-radiating and mainly confined within the jet core, featuring the larger wavenumber, and the other two consisting of toroidal and flapping Kelvin-Helmholtz signatures of identical wavelength and their associated downstream acoustic radiation, similarly to the modes discussed in the previous section.

The separation into $k_x < 0$ and $k_x > 0$ components, crucial in this case to understand the structure of the modes, is shown in figures 11 and 12 for the first and second SPOD modes, respectively. The splitting reveals a positive phase-speed component consisting of mixed $m = 0$ and $m = 1$ wavepacket structures, together with a negative phase-speed component formed by an $m = 1$ trapped acoustic wave with support within the supersonic jet-core region. Such acoustic wave, which strongly resembles
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Contours of the real part of the leading symmetric SPOD mode for $M_j = 1.26$ at $St = 0.56$: (left) $k_x < 0$ component; (right) $k_x > 0$ component; (top) $\sigma'$ field, $C_{\sigma'}$ SPOD problem; (middle) $\Theta'$ field, $C_{\Theta'}$ SPOD; (bottom) $\sigma'$ field, $C_{\Theta'}$ SPOD.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10** Contours of the real part of the antisymmetric SPOD modes for $M_j = 1.26$ at $St = 0.36$: (top) $\sigma'$ field, $C_{\sigma'}$ SPOD problem; (middle) $\Theta'$ field, $C_{\Theta'}$ SPOD; (bottom) $\sigma'$ field, $C_{\Theta'}$ SPOD; (left) mode 1; (right) mode 2.</figcaption></figure>

the so-called duct modes [67,19], is in this case an example of a downstream-propagating wave with negative phase velocity. Its positive group velocity can be estimated by performing the DFT along
<figure><img src="image_12.png" /><figcaption>**Fig. 11** Contours of the real part of the antisymmetric SPOD mode 1 for $M_j = 1.26$ at $St = 0.36$: (left) $k_x < 0$ component; (right) $k_x > 0$ component; (top) $\sigma'$ field, $C_{\sigma'}$ SPOD problem; (middle) $\Theta'$ field, $C_{\Theta'}$ SPOD; (bottom) $\sigma'$ field, $C_{\Theta'}$ SPOD.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12** Contours of the real part of the antisymmetric SPOD mode 2 for $M_j = 1.26$ at $St = 0.36$: (left) $k_x < 0$ component; (right) $k_x > 0$ component; (top) $\sigma'$ field, $C_{\sigma'}$ SPOD problem; (middle) $\Theta'$ field, $C_{\Theta'}$ SPOD; (bottom) $\sigma'$ field, $C_{\Theta'}$ SPOD.</figcaption></figure>

$x$ and tracking the evolution of its dominant streamwise wavenumber as a function of frequency, as illustrated in the $St-k_x$ maps found in [41].

The second and third rows of figures 10 to 12 show the performance of $C_{\Theta'}$ SPOD in separating the different coherent structures found at this frequency. The first $C_{\Theta'}$ SPOD mode comprises the toroidal wavepacket structure, together with a weak signature of the trapped acoustic wave. In turn, the second
$C_{\Theta'}$ SPOD mode consists of the flapping wavepacket together with a much stronger signature of the $m=1$ acoustic wave. These observations are further reflected on the reconstructed schlieren fields of the $C_{\Theta'}$ SPOD modes. Whereas the $C_{\sigma'}$ SPOD mixes all three $m=0$ and $m=1$ waves in the first and second modes, the $C_{\Theta'}$ SPOD tends to separate them effectively into different SPOD modes. The $m=0$ wave is recovered in the first SPOD mode, while the two $m=1$ waves are mainly found in the second SPOD mode, and are distinguished by the phase-speed splitting. Although the separation is not perfect, the latter behavior is more advantageous for the visualization and physical interpretation of the different type of oscillations undergone by the jets.

## 5 Conclusions

A methodology to extract coherent information from high-speed schlieren images of twin jets with enhanced interpretability has been presented. Following Prasad & Gaitonde [44], the proposed approach makes use of the momentum decomposition introduced by Doak [15], originally applied to a complex model problem by Jordan et al. [28], to derive a Poisson equation relating the schlieren fluctuation field with the streamwise gradient of the momentum potential fluctuations integrated along the line of sight, hereby denoted by $\Theta'$. As opposed to schlieren images, which are dominated in the turbulent mixing region by vortical fluctuations of a broad range of length scales, $\Theta'$ provides a representation of the potential (acoustic and thermal) energy embedded in the fluctuations, whose spatial structure in the frequency domain is very organized. By integrating $\Theta'$ within the spectral proper orthogonal decomposition (SPOD) framework, coherent structures associated with the dynamics of the momentum potential fluctuation field can be obtained. With respect to the technique of Prasad & Gaitonde [44], three aspects of the proposed methodology are novel:

- The Poisson equation is solved in the spectral (frequency-wavenumber) domain by means of the discrete Fourier transform, instead of the space-time domain. This simplifies the calculation of an accurate time derivative of the schlieren fluctuations, leveraging the necessity to transform the schlieren data to the frequency domain for the purpose of performing SPOD, and reduces the computational cost of solving the Poisson equation by exchanging the solution of a linear system of equations by fast Fourier transforms in each spatial direction.

- The impact of solving the Poisson equation in a truncated domain with *ad hoc* boundary conditions is assessed, illustrating that, in practice, unphysical harmonic waves are artificially introduced in the domain which contaminate the solution. A procedure is proposed to remove these spurious components in the spectral domain.

- The solution of the SPOD eigenvalue problem based on the cross-spectral density of $\Theta'$ is used to reconstruct schlieren SPOD modes, therefore obtaining coherent structures of the density gradient fluctuation field that obey the low-order dynamics of the momentum potential field instead of those of the schlieren field.

The methodology has been applied to schlieren visualizations of supersonic twin jets at perfectly-expanded and overexpanded conditions. When the SPOD is performed using $\Theta'$ realizations, coherent structures of the twin-jet field are obtained which follow a lower-rank classification than those obtained when the SPOD problem is based on schlieren images alone. The SPOD modes calculated when $\Theta'$ is used to construct the cross-spectral density tensor retain the coherent structure and the acoustic radiation associated with Kelvin-Helmholtz wavepackets and the screech resonance mechanism, but are substantially more effective in separating the toroidal ($m=0$) and flapping ($m=1$) fluctuations coexisting in coupled twin jets, which are recovered as independent SPOD modes. In addition, as the modal decomposition is based on a variable which represents potential-energy fluctuations, the extracted coherent structures are more closely related to the sound-generating mechanisms of the system and facilitate a more detailed analysis of the underlying physics. The comparison between $C_{\sigma'}$ SPOD modes and reconstructed $C_{\Theta'}$ schlieren SPOD modes demonstrates the higher interpretability of the latter structures in favor of physical understanding of jet noise mechanisms.

**Acknowledgements** The authors acknowledge S. Girard and D. Eysseric for the setup, technical assistance and critical assessment of the schlieren measurements at the PROMÉTÉE platform of Institut Pprime.
## Declarations

### Ethical Approval

Not applicable.

### Funding

The work of D.R. and I.P received partial financial support from the Government of the Community of Madrid within the multi-annual agreement with Universidad Politécnica de Madrid through the Program of Excellence in Faculty (V-PRICIT line 3) and the Program of Impulse of Young Researchers (V-PRICIT lines 1 and 3, Grant No. APOYO-JOVENES-WYOWRI-135-DZBLJU). I.P also received partial funding from the European Union's NextGenerationEU program (RD 289/2021: Ayudas Margarita Salas para la formación de jóvenes doctores).

## Availability of data and materials

The datasets generated and/or analyzed in this investigation are available from the corresponding author upon reasonable request.

## Appendix A Unphysical harmonic waves in $\Theta'$ and the impact of the filtering process used to remove them

As discussed in §3.4, the solution of the Poisson equation (8) to compute $\Theta'$ fields from schlieren fluctuations would require impractically large domains to reach boundaries at which the flow is irrotational, or appropriate boundary conditions for $\Theta'$. Due to the fact that, in general, it is not possible to determine the boundary values that would provide an irrotational potential field at the boundaries of the domain, an unphysical harmonic component is introduced in the solution when other boundary conditions (such as homogeneous Dirichlet or periodic conditions) are enforced. The solution of the Poisson problem in the spectral domain allows the identification of such harmonic component as spurious waves in the solution, and to filter them out based on their unphysical phase velocities.

To illustrate the impact of the unphysical harmonic waves on the calculated SPOD modes when they are not removed, figure 13 depicts the amplitude function of the first symmetric SPOD mode obtained for $M_j = 1.26$ at the fundamental screech frequency ($St = 0.56$). The left column shows the schlieren SPOD mode obtained from the schlieren-based SPOD problem, originally presented in figure 8 but repeated here for convenience, while the right column shows the $\Theta'$ SPOD mode obtained without filtering the unphysical energetic component. In contrast to the $\Theta'$ amplitude function shown in figure 8, when the energetic harmonic waves are not filtered from $\tilde{\Theta}'$, the $\Theta'$ SPOD mode is contaminated with small wavenumber (large wavelength) unphysical structures, as shown in the top right contour plot of figure 13. The energetic signature of these waves can be easily identified by performing the Fourier transform of the SPOD mode along $x$, as shown in the bottom right plot of the figure. For comparison, the energy distribution of the original schlieren mode is also displayed in the bottom left picture. As expected, the energy signature of this mode mainly consists of the Kelvin-Helmholtz energy band, located at positive $k_x$ slightly larger than the acoustic wavenumber, and an acoustic energy band resulting from the upstream- and downstream-traveling acoustic waves excited in the screech resonance, mainly located at negative wavenumbers below the acoustic one.

By comparison with the energy distribution of the original schlieren mode, the unfiltered $\Theta'$ mode contains an additional strong energetic signature centered at $k_x = 0$, which is not visible in the schlieren map. This energetic band is purely due to the artificial harmonic waves introduced in the numerical solution of the Poisson equation. Note that its signature can be distinguished from the physical energy components described above, and that it is composed of structures that travel at unphysical, high-supersonic speeds. This result explains the origin of the filtering procedure described in §3.4.

By removing the energetic components with wavenumbers below a chosen phase-speed threshold, in this case $c_{ph,c} = 1.2c_\infty$, the screech SPOD modes based on $C_{\Theta'}$ presented in figure 8 are obtained.
<figure><img src="image_14.png" /><figcaption>**Fig. 13** Symmetric SPOD mode 1 for $M_j = 1.26$ at $St = 0.56$: (left) results from the schlieren-based SPOD problem; (right) results from the $\Theta'$-based SPOD problem without filtering unphysical harmonic waves; (top) contours of the real part of the mode amplitude function; (bottom) Fourier transform of the SPOD mode along $x$, represented as contours of the logarithm of the energy magnitude. The red dash-dot lines denote the acoustic wavenumber associated to the frequency under investigation, given by the value at which the phase speed equals the freestream speed of sound. The white dashes lines represent the nozzle lip lines.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14** Contours of the streamwise Fourier transform of the first symmetric SPOD mode for $M_j = 1.26$ at $St = 0.56$, obtained from the $\Theta'$-based SPOD problem with filtered unphysical harmonic waves (see figure 8): (left) logarithm of the energy magnitude of the $\Theta'$ SPOD mode; (right) logarithm of the energy magnitude of the schlieren mode reconstructed from the $\Theta'$-based SPOD solution.</figcaption></figure>

The energy signature of those modes is shown in figure 14. Note that, now, the energy of the $\Theta'$ mode only shows the signatures of the physical waves playing a role in the screech resonance. In addition, the reconstructed schlieren mode has an energy distribution that is practically identical to the original schlieren SPOD mode, which demonstrates the effectiveness of the filtering approach and that it does not impact the energy of the physical information of interest.

Another case is also reported to justify the need to identify and control the unphysical harmonic waves introduced by the *ad hoc* imposition of boundary conditions in the Poisson equation. Figure 15 shows the unfiltered antisymmetric $\Theta'$ SPOD results for $M_j = 1.54$ at $St = 0.61$, in direct compari-
<figure><img src="image_16.png" /><figcaption>**Fig. 15** Antisymmetric SPOD modes for $M_j = 1.54$ at $St = 0.61$ obtained from the $\Theta'$-based SPOD calculation without filtering unphysical harmonic waves: (left) first SPOD mode; (right) second SPOD mode; (top) contours of the real part of the mode amplitude function; (bottom) Fourier transform of the SPOD mode along $x$, represented as contours of the logarithm of the energy magnitude.</figcaption></figure>

son with the filtered results already presented in figure 6. Both the first and second SPOD modes are depicted in the top row, once again showing the impact of the unphysical waves on the amplitude functions, which is more significant for the first mode in this case. This is corroborated by the distribution of the energy magnitude over the streamwise wavenumbers, illustrated in the bottom row. Two main energy bands are visible for each mode, the one associated with the Kelvin-Helmholtz waves, located very close to the acoustic wavenumber, and the one associated with the harmonic unphysical energy components, which is stronger for the first mode. Note that in this case no physical energy components are present at negative wavenumbers, which reflects the fact that the mixing noise associated with Kelvin-Helmholtz instabilities is the dominant mechanism in the perfectly-expanded regime.

The corresponding energy signatures for the filtered modes are shown in figure 16 for comparison, including the schlieren modes from the $C_{\sigma'}$ SPOD and the $C_{\Theta'}$ SPOD results. The removal of the unphysical harmonic waves from the $C_{\Theta'}$ SPOD problem yields modes whose energetic structure is consistent with that of the original schlieren modes, providing evidence that the applied filtering procedure is suitable for the problem under study.

## A.1 Effect of the phase-speed threshold employed for filtering the most energetic unphysical harmonic waves

The choice of the phase-speed threshold employed to remove the most energetic harmonic waves can have an impact on the SPOD eigenvalue problem. To assess the sensitivity of the SPOD results with respect to $c_{ph,c}$, three different threshold values are considered, namely $c_{ph,c} = 1.2c_{\infty}$, $1.5c_{\infty}$ and $2c_{\infty}$.

Figure 17 presents a comparison of the antisymmetric SPOD spectra ($M_j = 1.26$) for the first and second SPOD modes obtained using the different filtering thresholds. The unfiltered SPOD spectra are also included for comparison. Two different important observations can be made. On the one hand, the comparison between filtered and unfiltered results shows that the unphysical energy components have a stronger influence at higher frequencies. On the other hand, the filtered results are not very sensitive to the selected phase-speed threshold. Their sensitivity is higher for lower frequencies because
<figure><img src="image_17.png" /><figcaption>**Fig. 16** Contours of the streamwise Fourier transform of the symmetric SPOD modes for $M_j = 1.54$ at $St = 0.61$, presented as the logarithm of the energy magnitude: (left) first SPOD mode; (right) second SPOD mode; (top) schlieren SPOD mode obtained from the schlieren-based SPOD problem; (middle) $\Theta'$ SPOD mode obtained from the $\Theta'$-based SPOD problem with filtered unphysical harmonic waves; (bottom) schlieren mode reconstructed from the $\Theta'$-based SPOD solution with filtered unphysical harmonic waves. The contour plots in this figure have a one-to-one correspondence with the amplitude functions displayed in figure 6.</figcaption></figure>

the acoustic wavenumber decreases with frequency. For very low frequencies, the unphysical harmonic energy band and the physical energy bands become very close together. Therefore, small values of $c_{ph,c}$ can impact the physical energy bands as well, which explains the small effect shown in the spectra. In particular, at low frequencies, small phase-speed thresholds can have a non-negligible effect on the growth of the coherent structures (wavepackets), which is governed by small wavenumbers dictated by the wavepacket envelope shape. For high frequencies, however, the range of wavenumbers contaminated by the artificial harmonic components is far from the acoustic wavenumber where the physical energy is generally concentrated. As a result, the spectra are very insensitive to the filter width in this regime.

For the results presented in this work, the influence of the employed phase-speed threshold ($1.2c_\infty$) is negligible. Figures 18 and 19 illustrate this by respectively comparing the first and second SPOD modes at $St = 0.36$ (the lowest frequency analyzed in the results section) for the three different filtering thresholds under consideration. The impact of the filter width on the amplitude function and the corresponding energy signature is very small; only the growth rate of the coherent structures of the second SPOD mode shows small deviations with different filter width. For the analysis of lower-frequency coherent structures, however, larger values of $c_{ph,c}$ are recommended.
<figure><img src="image_18.png" /><figcaption>**Fig. 17** Comparison of antisymmetric $\Theta'$ SPOD spectra for $M_j = 1.26$ using different phase-speed thresholds for removing unphysical harmonic waves: (left) first SPOD mode; (right) second SPOD mode. The unfiltered case is also added for comparison.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 18** Comparison of the first antisymmetric $\Theta'$ SPOD mode for $M_j = 1.26$, $St = 0.36$ using different phase-speed thresholds for removing unphysical harmonic waves: (left) $c_{ph,c} = 1.2c_\infty$; (middle) $c_{ph,c} = 1.5c_\infty$; (right) $c_{ph,c} = 2c_\infty$; (top) real part of the SPOD mode amplitude function; (bottom) logarithm of the energy magnitude from the streamwise Fourier transform of the SPOD mode. The light blue dash-dot lines denote the wavenumber associated with the chosen phase-speed threshold for each case.</figcaption></figure>

## Appendix B Comparison of alternative approaches to construct the $\Theta'$ SPOD eigenvalue problem

Different approaches can be considered to construct the SPOD eigenvalue problem based on the cross-spectral density of $\Theta'$. Besides the methodology outlined in section 3.4, where the Poisson equation is solved in the frequency-wavenumber domain, two other procedures have also been tested, which are summarized below.
<figure><img src="image_20.png" /><figcaption>**Fig. 19** Comparison of the second antisymmetric $\Theta'$ SPOD mode for $M_j = 1.26$, $St = 0.36$ using different phase-speed thresholds for removing unphysical harmonic waves: (left) $c_{ph,c} = 1.2c_\infty$; (middle) $c_{ph,c} = 1.5c_\infty$; (right) $c_{ph,c} = 2c_\infty$; (top) real part of the SPOD mode amplitude function; (bottom) logarithm of the energy magnitude from the streamwise Fourier transform of the SPOD mode.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 20** Comparison of the symmetric $\Theta'$ SPOD spectra for $M_j = 1.26$ obtained with different approaches: (top-left) using $\Theta'$ snapshots computed in the space-time domain; (top-right) using $\Theta'$ computed in the space-frequency domain; (bottom-left) using $\Theta'$ computed in the wavenumber-frequency domain; (bottom-right) comparison of the sum of the spectral density of all the SPOD modes for the different approaches.</figcaption></figure>

## B.1 Calculation of $\Theta'$ snapshots in the space-time domain

In this approach, $\Theta'$ snapshots are directly computed in the space-time domain through the numerical solution of equation (8). The time derivative of schlieren fluctuations (source term of the Poisson
equation) is in this case computed by means of a central finite difference scheme, for which different orders of accuracy have been tested. Similarly, the Laplacian operator is discretized using a second-order central finite-difference scheme for a uniform grid. Regarding the boundary conditions, an homogeneous Dirichlet condition is imposed. This results in a linear system of equations which is solved by means of sparse Cholesky factorization using the CHOLMOD package [9]. The SPOD algorithm is then applied directly on $\Theta'$ snapshots. Once the realizations of $\Theta'$ are Fourier-transformed in time, the FFT of $\tilde{\Theta}'$ along $x$ is calculated and the most energetic harmonic components are removed in the same way as introduced before.

The symmetric SPOD spectrum resulting from this approach is represented in the top left plot of figure 20 for $M_j = 1.26$, obtained using a second-order scheme for computing the time derivative of the schlieren fluctuations. Although the broadband mixing noise region and the screech resonance tones are well predicted, this method leads to a very fast decay of the SPOD energy for high frequencies. The reason for this lies on the low-order of the scheme employed for the calculation of the time derivative of the schlieren fluctuations, which attenuates the high-frequency content present in the original snapshots. In order to preserve the energy at high frequencies, the time derivative should be computed by means of a high-order discretization scheme in the time domain. To illustrate this, the total spectral density obtained with fourth and eighth-order schemes for the time derivative are also shown in the bottom right plot of figure 20. As expected, higher-order time derivative calculations preserve a larger portion of the high frequency content, although an order significantly higher than 8 would be required to correctly approximate the energy of the highest frequencies. In the limit in which a spectral discretization scheme (order equal to the number of instantaneous snapshots) is considered, the time-domain calculation of $\Theta'$ would theoretically recover the same SPOD spectrum as the frequency-domain approach.

In practice, standard central finite-difference schemes have accuracy limitations when the order becomes very high for a given point spacing, owing to the fact that the stencils become ill-conditioned [52]. To achieve very high-order differentiation in the time domain, the use of more elaborate discretization techniques is required, such as Padé finite difference schemes, complex-variable differentiation or Chebyshev spectral collocation methods [37,70]. These techniques, however, increase the computational complexity and the implementation complexity of the algorithm. For this reason, the calculation of the time derivative in the frequency domain is argued to be more convenient for the methodology presented in this study.

The application of homogeneous Dirichlet boundary conditions in this approach may be combined with the use of a sponge layer to force the schlieren fluctuations to decay to zero towards the boundaries and avoid discontinuities (as done in [44,45]). Here, the use of a Planck-taper spatial window has been tested for that purpose. In all the cases investigated, using tapering parameter values of $\varepsilon = 0.1, 0.15$ and $0.2$, no significant reduction of the unphysical harmonic energetic components introduced in the solution could be achieved. Therefore, the filtering process is still required to obtain meaningful coherent structures.

As shown in Appendix C, the computational cost of this approach is higher than the SPOD methodology presented in this work, which computes $\Theta'$ realizations in the wavenumber-frequency domain.

## B.2 Calculation of $\tilde{\Theta}'$ in the space-frequency domain

In this case, snapshots of $\Theta'$ are not obtained in the time domain. Instead, the calculation of $\tilde{\Theta}'$ is carried out by solving the Poisson equation in the space-frequency domain. This is achieved by computing the time derivative of schlieren fluctuations in the frequency domain, but solving the resulting Poisson equation (11) in the $x, y$ spatial domain using a second-order finite-difference scheme, as in the previous approach.

The spectrum obtained via this method is displayed in the top right plot of figure 20. With this approach, the obtained spectrum is almost identical to the one obtained solving for $\Theta'$ in the wavenumber-frequency domain. However, note that this procedure also requires the solution of a linear system of equations for each $\Theta'$ field and its posterior spatial Fourier transform along $x$ to remove the unphysical harmonic waves. Therefore, the computational cost of this approach is very similar to the space-time domain variant presented in the previous section. Its advantage over the previous one is the recovery of the high frequency content thanks to the spectral calculation of the time derivative.
<figure><img src="image_22.png" /><figcaption>**Fig. 21** Comparison of the computation time $\tau$ between different approaches as a function of the number of snapshots: (left) total computation time in minutes; (right) computation time of the two $C_{\Theta'}$ methods relative to the computation time of $C_{\sigma'}$.</figcaption></figure>

## Appendix C Computational cost comparison between the space-time and the wavenumber-frequency domain approaches for the $\Theta'$ SPOD problem

To provide a quantitative comparison between the computational cost of the space-time $(x, y, t)$ and the wavenumber-frequency $(k_x, k_y, \omega)$ domain approaches for solving the SPOD problem based on $\Theta'$, figure 21 reports the computation time $(\tau)$ resulting from both approaches as a function of the number of snapshots. The computation time of the schlieren SPOD problem ($C_{\sigma'}$) is also added for reference. The computation time is reduced by a factor of 2.5 when the Poisson equation is solved in the spectral domain instead of the space-time domain. With respect to the SPOD based on schlieren snapshots, the $C_{\Theta'}(k_x, k_y, \omega)$ approach increases the computation time by 1.5, while the $C_{\Theta'}(x, y, t)$ method increases it by 3.8.

The SPOD calculations have been performed using MATLAB R2021b on a 64-bit Linux system with a AMD Ryzen 5 5600X CPU. The MATLAB built-in functions for computing the FFT are used, which employ the FFTW [22] library. Similarly, the built-in functions for solving the eigenvalue problem are also employed, which use LAPACK [1] routines. In the space-time approach, the linear system of equations resulting from the spatial discretization of the Poisson equation is solved by means of sparse Cholesky factorization via the CHOLMOD [9] library. For the studied cases, this solution method was found to be faster than the different MATLAB built-in implementations of the conjugate gradient iterative method.

The following aspects have also been considered for the comparison:

- In the space-time $\Theta'$-SPOD algorithm, the overlap between blocks is taken into account to avoid computing the same $\Theta'$ realizations more than one time.

- In the space-time $\Theta'$-SPOD algorithm, second-order finite difference discretizations are employed both in space and time.

- In the wavenumber-frequency $\Theta'$-SPOD algorithm, the original schlieren domain (352 x 512 points) is zero-padded to 1024 x 1024 points.

The computational performance reported in this appendix is subject to change depending on the implementation details and optimization of each algorithm. However, the significant speedup obtained in the $(k_x, k_y, \omega)$ approach justifies the advantage of solving the Poisson equation in the spectral domain for SPOD purposes.

## References

1. Anderson, E. Bai, Z. Bischof, C. Blackford, S. Demmel, J. Dongarra, J. Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., Sorensen, D.: LAPACK Users' Guide, third edn. Society for Industrial
and Applied Mathematics, Philadelphia, PA (1999)

2. Batchelor, G.K.: An Introduction to Fluid Dynamics. Cambridge Mathematical Library. Cambridge University Press, Cambridge (2000)

3. Berkooz, G., Holmes, P., Lumley, J.L.: The proper orthogonal decomposition in the analysis of turbulent flows. Annu. Rev. Fluid Mech. **25**(1), 539–575 (1993)

4. Bhat, W.V.: Acoustic characteristics of two parallel flow jets. AIAA Paper 77-1290 (1977)

5. Borée, J.: Extended proper orthogonal decomposition: A tool to analyse correlated events in turbulent flows. Experiments in Fluids **35**(2), 188–192 (2003)

6. Cavalieri, A., Jordan, P., Colonius, T., Gervais, Y.: Axisymmetric superdirectivity in subsonic jets. J. Fluid Mech. **704**, 388–420 (2012)

7. Cavalieri, A., Rodríguez, D., Jordan, P., Colonius, T., Gervais, Y.: Wavepackets in the velocity field of turbulent jets. J. Fluid Mech. **730**, 559–592 (2013)

8. Cavalieri, A.V.G., Daviller, G., Comte, P., Jordan, P., Tadmor, G., Gervais, Y.: Using large eddy simulation to explore sound-source mechanisms in jets. J. Sound Vib. (330), 4098–4113 (2011)

9. Chen, Y., Davis, T.A., Hager, W.W., Rajamanickam, S.: Algorithm 887: Cholmod, supernodal sparse cholesky factorization and update/downdate. ACM Trans. Math. Softw. **35**(3) (2008)

10. Chu, B.T.: On the energy transfer to small disturbances in fluid flow (Part I). Acta Mechanica **1**(3), 215–234 (1965)

11. Citriniti, J.H., George, W.K.: Reconstruction of the global velocity field in the axisymmetric mixing layer utilizing the proper orthogonal decomposition. J. Fluid Mech. **418**, 137–166 (2000)

12. Crighton, D.G., Gaster, M.: Stability of slowly diverging jet flow. J. Fluid Mech. **77**(2), 397–413 (1976)

13. Crighton, D.G., Huerre, P.: Shear-layer pressure fluctuations and superdirective acoustic sources. J. Fluid Mech. **220**, 355–368 (1990)

14. Crow, S., Champagne, F.: Orderly structure in jet turbulence. J. Fluid Mech. **48**(3), 547–591 (1971)

15. Doak, P.E.: Momentum potential theory of energy flux carried by momentum fluctuations. Journal of Sound and Vibration **131**(1), 67–90 (1989)

16. Du, Z.: Acoustic and kelvin-helmholtz instability waves of twin supersonic jets. Ph.D. thesis, The Florida State University (1993)

17. Edgington-Mitchell, D., Jaunet, V., Jordan, P., Towne, A., Soria, J., Honnery, D.: Upstream-travelling acoustic jet modes as a closure mechanism for screech. J. Fluid Mech. **855**, 1–12 (2018)

18. Edgington-Mitchell, D., Li, X., Liu, N., He, F., Wong, T.Y., MacKenzie, J., Nogueira, P.: A unifying theory of jet screech. J. Fluid Mech. **945**, 1–24 (2022)

19. Edgington-Mitchell, D., Wang, T., Nogueira, P., Schmidt, O., Jaunet, V., Duke, D., Jordan, P., Towne, A.: Waves in screeching jets. Journal of Fluid Mechanics **913** (2021)

20. Esfahani, A., Webb, N.J., Samimy, M.: Control of coupling in twin rectangular supersonic jets. AIAA Aviation Forum 2021, Virtual event, AIAA Paper 2021-2122 (2021)

21. Freund, J.B., Colonius, T.: POD analysis of sound generation by a turbulent jet. 40th AIAA Aerospace Sciences Meeting and Exhibit, January 14-17, 2002, Reno, NV, AIAA Paper 2002-0072 (2002)

22. Frigo, M., Johnson, S.G.: The design and implementation of FFTW3. Proceedings of the IEEE **93**(2), 216–231 (2005). Special issue on “Program Generation, Optimization, and Platform Adaptation”

23. Garnaud, X., Lesshafft, L., Schmid, P.J., Huerre, P.: The preferred mode of incompressible jets: linear frequency response analysis. J. Fluid Mech. **716**, 189–202 (2013)

24. Goparaju, K., Gaitonde, D.V.: Dynamics of Closely Spaced Supersonic Jets. Journal of Propulsion and Power **34**(2), 327–339 (2018)

25. Gudmundsson, K., Colonius, T.: Instability wave models for the near-field fluctuations of turbulent jets. J. Fluid Mech. **689**, 97–128 (2011)

26. Jeun, J., Nichols, J.W., Jovanović, M.R.: Input–output analysis of high-speed axisymmetric isothermal jet noise. Phys. Fluids **28**(4), 047,101 (2016)

27. Jordan, P., Colonius, T.: Wave Packets and Turbulent Jet Noise. Annu. Rev. Fluid Mech. **45**(1), 173–195 (2013)

28. Jordan, P., Daviller, G., Comte, P.: Doak’s momentum potential theory of energy flux used to study a solenoidal wavepacket. Journal of Sound and Vibration **332**(17), 3924–3936 (2013)

29. Juvé, D., Sunyach, M., Compte-Bellot, G.: Intermittency in the noise emission in subsonic cold jets. J. Sound Vib. **71**, 319–332 (1980)

30. Kantola, R.A.: Acoustic properties of heated twin jets. J. Sound and Vibration **79**(1), 79–106 (1981)

31. Kaplan, O., Jordan, P., Cavalieri, A.V., Brès, G.A.: Nozzle dynamics and wavepackets in turbulent jets. Journal of Fluid Mechanics **923**, 1–27 (2021)

32. Karban, U., Bugeat, B., Towne, A., Lesshafft, L., Agarwal, A., Jordan, P.: An empirical model of noise sources in subsonic jets. Journal of Fluid Mechanics **965** (2023)

33. Karban, U., Martini, E., Cavalieri, A.V., Lesshafft, L., Jordan, P.: Self-similar mechanisms in wall turbulence studied using resolvent analysis. Journal of Fluid Mechanics **939**, 1–31 (2022)

34. Karnam, A., Saleem, M., Gutmark, E.: Influence of nozzle geometry on screech instability closure. Physics of Fluids **35**(8) (2023)

35. McKechan, D.J., Robinson, C., Sathyaprakash, B.S.: A tapering window for time-domain templates and simulated signals in the detection of gravitational waves from coalescing compact binaries. Classical and Quantum Gravity **27**(8) (2010)

36. Michalke, A.: Survey on jet instability theory. Prog. Aero. Sci. **21**, 159–199 (1984)

37. Moin, P.: Fundamentals of engineering numerical analysis, second edn. Cambridge University Press, New York (2010)

38. Morris, P.J.: Instability waves in twin supersonic jets. J. Fluid Mech. **220**, 293–307 (1990)
39. Nogueira, P., Stavropoulos, M., Edgington-Mitchell, D.: Wavepacket coupling in screeching twin-jets. Annual Conference of the Australian Acoustical Society 2021 pp. 123–130 (2021)

40. Nogueira, P.A., Edgington-Mitchell, D.M.: Investigation of supersonic twin-jet coupling using spatial linear stability analysis. Journal of Fluid Mechanics **918**, A38 (2021)

41. Padilla-Montero, I., Rodríguez, D., Jaunet, V., Girard, S., Eysseric, D., Jordan, P.: Investigation of coherent motions and noise radiation in twin supersonic jets using high-speed Schlieren images. 10th Convention of the European Acoustics Association, September 11-15, 2023, Turin, Italy (2023)

42. Padilla-Montero, I., Rodríguez, D., Jaunet, V., Girard, S., Eysseric, D., Jordan, P.: Mean flow and linear stability analysis of twin supersonic jets. AIAA Aviation Forum 2023, June 12-16, 2023, San Diego, CA, AIAA Paper 2023-3350 (2023)

43. Piot, E., Casalis, G., Muller, F., Bailly, C.: Investigation of the PSE approach for subsonic and supersonic hot jets. Detailed comparisons with LES and Linearized Euler Equations results. Int. J. Aeroacoustics **5**, 361–393 (2006)

44. Prasad, C., Gaitonde, D.V.: A robust physics-based method to filter coherent wavepackets from high-speed schlieren images. Journal of Fluid Mechanics **940**, 1–11 (2022)

45. Prasad, C., Gaitonde, D.V., Esfahani, A., Webb, N.J., Samimy, M.: Examination of wavepackets in forced and unforced rectangular twin jets with high-speed schlieren. AIAA Scitech Forum 2022, January 3-7, 2022, San Diego, CA, AIAA Paper 2022-2402 (2022)

46. Prasad, C., Morris, P.J.: Steady active control of noise radiation from highly heated supersonic jets. The Journal of the Acoustical Society of America **149**(2), 1306–1317 (2021)

47. Ray, P., Cheung, L., Lele, S.: On the growth and propagation of linear instability waves in compressible turbulent jets. Physics of Fluids **21**, 054,106 (2009)

48. Rodríguez, D.: Wavepacket models for supersonic twin-jets. In: AIAA AVIATION 2021 FORUM, pp. 1–16. American Institute of Aeronautics and Astronautics, Reston, Virginia (2021)

49. Rodríguez, D., Jotkar, M.R., Gennaro, E.M.: Wavepacket models for subsonic twin jets using 3D parabolized stability equations. Comptes Rendus Mécanique **346**(10), 890–902 (2018)

50. Rodríguez, D., Sinha, A., Brès, G.A., Colonius, T.: Inlet conditions for wave packet models in turbulent jets based on eigenmode decomposition of large eddy simulation data. Physics of Fluids **25**(10), 105,107 (2013)

51. Rodríguez, D., Stavropoulos, M.N., Nogueira, P.A., Edgington-Mitchell, D.M., Jordan, P.: On the preferred flapping motion of round twin jets. Journal of Fluid Mechanics **977**, A4 (2023)

52. Sauer, T.: Numerical Analysis, second edn. Pearson Education, Boston (2012)

53. Schmidt, O.T., Colonius, T.: Guide to spectral proper orthogonal decomposition. AIAA Journal **58**(3), 1023–1033 (2020)

54. Schmidt, O.T., Towne, A., Rigas, G., Colonius, T., Brès, G.: Spectral analysis of jet turbulence. J. Fluid Mech. **855**, 953–982 (2018)

55. Schoder, S., Roppert, K., Kaltenbacher, M.: Helmholtz’s decomposition applied to aeroacoustics. 25th AIAA/CEAS Aeroacoustics Conference, 2019 (May), 1–11 (2019)

56. Schoder, S., Roppert, K., Kaltenbacher, M.: Helmholtz’s decomposition for compressible flows and its application to computational aeroacoustics. Partial Differential Equations and Applications **1**(6), 1–20 (2020)

57. Sedel'nikov, T.K.: The dispersion relations for multilayer jets and for several jets. Physics of Aerodynamic Noise (ed. A. V. Rimskiy-Korsakov, transl. NASA TTF-538) (1967)

58. Sinha, A., Rodríguez, D., Brès, G., Colonius, T.: Wavepacket models for supersonic jet noise. J. Fluid Mech. **742**, 71–95 (2014)

59. Sirovich, L.: Turbulence and the dynamics of coherent structures. part i: Coherent structures. Quarterly of Applied Mathematics **45**(3), 561–571 (1987)

60. Souza, D.S., Rodríguez, D., Himeno, F.H.T., Medeiros, M.: Dynamics of the large-scale structures and associated noise emission in airfoil slats. J. Fluid Mech. **875**, 1004–1034 (2019)

61. Souza, D.S., Rodríguez, D., Simões, L.G.C., Medeiros, M.A.F.: Effect of an excrescence in the slat cove: Flow-field, acoustic radiation and coherent structures. Aerospace Science and Technology **44**, 108–115 (2015)

62. Sprössig, W.: On helmholtz decompositions and their generalizations—an overview. Mathematical Methods in the Applied Sciences **33**(4), 374–383 (2010)

63. Stavropoulos, M.N., Mancinelli, M., Jordan, P., Jaunet, V., Weightman, J., Edgington-Mitchell, D.M., Nogueira, P.A.: The axisymmetric screech tones of round twin jets examined via linear stability theory. Journal of Fluid Mechanics **965**, 1–29 (2023)

64. Suzuki, T., Colonius, T.: Instability waves in a subsonic round jet detected using a near-field phased microphone array. J. Fluid Mech. **565**, 197–226 (2006)

65. Tam, C.K.W.: Supersonic jet noise. Annu. Rev. Fluid Mech. **27**(1), 17–43 (1995)

66. Tam, C.K.W., Hu, F.Q.: On the three families of instability waves of high-speed jets. Journal of Fluid Mechanics **201**, 447–483 (1989)

67. Towne, A., Cavalieri, A.V., Jordan, P., Colonius, T., Schmidt, O., Jaunet, V., Brès, G.A.: Acoustic resonance in the potential core of subsonic jets. Journal of Fluid Mechanics **825**, 1113–1152 (2017)

68. Towne, A., Colonius, T.: One-way spatial integration of hyperbolic equations. J. Comp. Phys. **300**, 844–861 (2015)

69. Towne, A., Schmidt, O.T., Colonius, T.: Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. Journal of Fluid Mechanics **847**, 821–867 (2018)

70. Trefethen, L.N.: Spectral methods in MATLAB. SIAM (2000)
71. Unnikrishnan, S., Gaitonde, D.V.: Acoustic, hydrodynamic and thermal modes in a supersonic cold jet. Journal of Fluid Mechanics **800**, 387–432 (2016)

72. Unnikrishnan, S., Gaitonde, D.V.: Transfer mechanisms from stochastic turbulence to organized acoustic radiation in a supersonic jet. European Journal of Mechanics, B/Fluids **72**, 38–56 (2018)

73. Wong, T.Y., Stavropoulos, M.N., Beekman, J.R., Towne, A., Nogueira, P.A., Weightman, J., Edgington-Mitchell, D.: Steady and unsteady coupling in twin weakly underexpanded round jets. Journal of Fluid Mechanics **964**, 1–43 (2023)