On the preferred flapping motion of round twin jets

Abstract Linear stability theory (LST) is often used to model the large-scale flow structures in the turbulent mixing region and near pressure field of high-speed jets. For perfectly expanded single round jets, these models predict the dominance of azimuthal wavenumbers $m=0$ and $m = 1$ helical modes for the lower frequency range, in agreement with empirical data. When LST is applied to twin-jet systems, four solution families appear following the odd/even behaviour of the pressure field about the symmetry planes. The interaction between the unsteady pressure fields of the two jets also results in their coupling. The individual modes of the different solution families no longer correspond to helical motions, but to flapping oscillations of the jet plumes. In the limit of large jet separations, when the jet coupling vanishes, the eigenvalues corresponding to the $m=1$ mode in each family are identical, and a linear combination of them recovers the helical motion. Conversely, as the jet separation decreases, the eigenvalues for the $m=1$ modes of each family diverge, thus favouring a particular flapping oscillation over the others and preventing the appearance of helical motions. The dominant mode of oscillation for a given jet Mach number $M_j$ and temperature ratio $T_R$ depends on the Strouhal number $St$ and jet separation $s$. Increasing both $M_j$ and $T_R$ independently is found to augment the jet coupling and modify the $(St,s)$ map of the preferred oscillation mode. Present results predict the preference of two modes when the jet interaction is relevant, namely varicose and especially sinuous flapping oscillations on the nozzles’ plane.


Introduction
Tactical fighters developed since the 1960s predominately feature fuselage-embedded twin jet engines.Additionally, multi-tube nozzles have been investigated as possible jet noise suppressors, and designs of future distributed-propulsion systems involve placing two or more parallel jet streams in close proximity.The closely spaced jets can interact both at the hydrodynamic and acoustic levels, giving rise to complex flow structures when compared with single round jets at equivalent operating conditions.
Following the seminal works of Mollo-Christensen (1967), the relation between the dominant components of the far-field noise radiated by high-speed jets and large-scale fluctuations in the turbulent mixing region, coherent over several nozzle diameters, has been the subject of growing research (Jordan & Colonius 2013;Cavalieri et al. 2019).Radiated sound is correlated with largescale, low-frequency fluctuations in the mixing region and with very few azimuthal modes for single isolated jets (Juvé et al. 1980;Hileman et al. 2005;Cavalieri et al. 2011).The existence of coherent structures in turbulent jets was first identified by Crow & Champagne (1971) and their resemblance to instability waves for harmonically-forced supersonic jets suggested the use of linear instability analysis to model them (Crighton & Gaster 1976;Michalke 1984).The presence of wavepackets in high-speed jets was finally demonstrated over the last two decades, as well as the ability of linear stability calculations to model them faithfully (Suzuki & Colonius 2006;Gudmundsson & Colonius 2011;Cavalieri et al. 2013).
As opposed to single round jets, instability analyses for twin jets are scarce in the literature, due to the mathematical complexity of the latter.The mean turbulent flow corresponding to isolated round jets is axially symmetric, enabling the introduction of azimuthal Fourier modes (each one characterised by an integer wavenumber ) and only requires spatial discretisation in the radial direction.Bipolar coordinates were used by Morris (1990) to study the inviscid instability of two axially-homogeneous parallel jets.He identified the counterparts of the different azimuthal Fourier modes known for single jets and classified them according to the symmetries about the jetcenter plane and the plane normal to it.Green & Crighton (1997) used a similar approach to study coupled oscillation modes of the jet cores, considering varicose and sinuous flapping motions of the two jets.More recently, Rodríguez et al. (2018) and Nogueira & Edgington-Mitchell (2021) analysed the local linear instability of twin-jet configurations by applying two-dimensional crossstream discretisations that do not restrict the spatial structure of the wavepackets.Interestingly, the latter analyses recover the same families of eigenmodes corresponding to the Fourier modes of single jets, albeit modified on account of the jet-jet interaction; additional eigenmodes corresponding to mechanisms not present for single jets were not identified.It was observed (but not reported) that the eigenfunctions for > 0 modes do not correspond to helical oscillations, but constitute combinations of the respective + and − helical modes with the precise phase relationship such that they describe flapping motions on the lateral or vertical planes.This is consistent with experimental observations and simulations in supersonic twin round jets.Seiner et al. (1988) and Wlezien (1989) showed that the B-mode of oscillation associated with screech manifests as a coupled flapping motion of the plumes occurring in the plane containing the jets.Alkislar et al. (2005) reported that the coupling between the two jets, for their particular conditions of jet Mach number and jet separation , results in a symmetric flapping with respect to the mid plane and dominant jet oscillations occuring in the plane containing the jets (referred to as lateral flapping).A similar observation was made by Kuo et al. (2017) based on Large Eddy Simulations.The experiments by Kuo et al. (2017) considered a twin-jet configuration with separation / = 2 and convergent-divergent nozzles with design Mach number of 1.23 and jet Mach numbers between = 1.15 and 1.4.The jets exhibited coupled oscillations, with a varicose (symmetric) lateral flapping motion being observed at conditions in which a single jet would present the helical oscillations typical of the screech B-modes.Experiments undertaken by Knast et al. (2018);Bell et al. (2018) reported both helical and flapping oscillations in twin jets, but due to experimental constraints the presence or absence of helical modes could not be rigorously determined.Later, Bell et al. (2021) demonstrated that, even at a fixed operation condition, twin-jet oscillations present intermittency.At some time lapses the motion of the two jets can be uncorrelated and present helical motions; at other lapses they are strongly correlated and present a coupled lateral flapping.
This paper revisits the locally-parallel linear instability of twin-jet configurations.The first objective is to show that the coupling of the two jets favours flapping motions over helical ones.The analysis is focused on = 1 modes, as experiments show their prevalence over other modes for most flow conditions in the supersonic regime, specially in the first diameters from the nozzle lip.The second objective is to map the preferred flapping mode for each Strouhal number and jet separation, and the impact of the jet Mach number and temperature ratio upon them.Two independent formulations of the linearised equations are used: (i) an approach that discretises the cross-stream plane using Cartesian coordinates, valid for of finite-thickness jets of arbitrary shape; and (ii) an inviscid vortex-sheet method analogous to that used by Morris (1990);Du (1993);Stavropoulos et al. (2021).Section 2 describes the two formulations.Section 3 presents the results of the analyses.The relevant eigenmode families are described briefly.The impact of the interaction of the pressure fields of the two jets on the eigenmodes and the appearance of preferred modes is then discussed.A parametric study is then presented that analyses the effect of the jet Mach number and temperature ratio on the preferred oscillation mode.The main conclusions are summarised in section 4.

Formulations of the linear stability problem for twin jets
The twin-jet geometry and geometrical parameters are shown in figure 1.The nozzle diameter is and the separation between the jet axes is .The streamwise coordinate is oriented perpendicular to the paper towards the reader.Radial and azimuthal coordinates measured from the origin are denoted by ( , ), and those measured with respect to the axis of each jet are denoted by ( 1 , 1 ) and ( 2 , 2 ).Physical quantities are made dimensionless using and the free-stream sound speed ∞ and density ∞ .Pressure is scaled with ∞ 2 ∞ and temperature with ( − 1) ∞ , where is the ratio of specific heats.The jet Mach number is defined as = / ∞ , with being the jet exit velocity.The jet temperature ratio is defined as = / ∞ , where is the jet exit temperature and ∞ the free-stream temperature.
Let ′ be a vector containing all the fluid variables of interest, e.g. the velocity = ( , , ), density , pressure and temperature .Linear instability theory (LST) studies small-amplitude disturbances superimposed on a time-invariant flow, either steady laminar or stationary turbulent mean flow, denoted here by ¯ .Invoking the locally-parallel flow assumption, modal disturbances of the form ′ ( , ) = ( , ) e i( − ) + c.c., (2.1) are introduced, where is the circular frequency, the streamwise wavenumber, the dimensionless time and c.c. denotes the complex conjugate.The Strouhal number defined as = / is related to the dimensionless angular frequency by = 2 . The derivation of the LST problem continues by substituting the modal decomposition (2.1) in the linearised governing equations and recasting the result as an eigenvalue problem.The spatial instability framework is used in this work, which consists of prescribing a real frequency and obtaining the corresponding eigensolutions as the eigenvalue/eigenfunction pairs ( , ).Thus, LST provides the dispersion relation D ( , ) = 0 that governs the evolution of linear instability waves.
LST is applied here to twin-jet configurations.As opposed to single round jets, axial symmetry cannot be exploited here to simplify the LST problem; however, the twin-jet mean flow is symmetric with respect to the ( , ) and ( , ) planes.Following Morris (1990), the LST solutions are separated in four families corresponding to the even or odd character of their pressure field with respect to the two planes.The two-letter notation used by Rodríguez et al. (2018); Nogueira & Edgington-Mitchell ( 2021) is adopted to classify the mode families, which is outlined in table 1.Two different formulations of the LST problem are used in this work, as described next.
T 1. Classification of mode families depending on the symmetries.The fourth column shows the relation of the notation used here with that by Morris (1990)

Formulation 1: Complete compressible Navier-Stokes equations discretised in cartesian coordinates
This formulation is valid for finite-thickness jets of arbitrary shape and has been applied previously to twin round jets (Rodríguez et al. 2018;Rodríguez 2021) and rectangular jets (Rodríguez et al. 2021).The viscous compressible continuity, momentum and energy equations in cartesian coordinates are used as the departure point.Ideal gas is assumed and the impact of temperature on viscosity is neglected.In the linearisation, density and temperature are used for the mean flow and pressure and temperature for the disturbances.The modal form (2.1) is introduced and the resulting generalised eigenvalue problem is recast in the form where matrix operators L and R depend on the mean flow and its derivatives, , the physical parameters , , and and parameters describing the twin-jet configuration (e.g./ , , ).The rectangular domain Ω = [0, ∞ ] × [0, ∞ ] is used in the discretisation.An analytical coordinate transformation is applied independently to the and coordinates to concentrate points around the jet shear layer.The odd/even behaviours of each family are imposed as boundary conditions along the = 0 and = 0 axes.For modes that are respectively symmetric (S) and anti-symmetric (A) with respect to the ( , )−plane, the conditions S : are imposed at = 0, and similarly for the behaviour with respect to the ( , )−plane.Vanishing of the disturbance velocity and temperature is imposed as far-field boundary conditions.A Neumann condition is imposed for the pressure.The rectangular domain Ω = [0, 7.5] × [0, 5] was checked to be large enough for the convergence of the results for the case with the largest jet separation considered herein ( / = 5), and is used for all calculations.
The linear operators L and R are discretised using a combination of variable-stencil highorder finite differences and sparse algebra that exploits the banded structure of the differentiation matrices.A 7-point stencil is used, which results in the optimal balance between accuracy and computational cost (Gennaro et al. 2013;Rodríguez & Gennaro 2017).After discretisation of the linear operators, the matrix eigenvalue problem (2.2) is solved using an in-house sparse implementation of the shift-and-invert Arnoldi's algorithm (Arnoldi 1951).Arnoldi's algorithm requires the solution of a number of linear problems, which is accomplished using the package MUMPS (Amestoy et al. 2001).
This formulation can also be applied without exploiting/imposing the symmetries, as done in Rodríguez et al. (2018).In this case, the computational domain used is and the same mode families are recovered.The results of the present formulation have been cross-validated with the formulation presented by Nogueira & Edgington-Mitchell (2021), that employs a Floquet formalism on the azimuthal direction to reduce the computation to a sector of the azimuthal domain and discretises it using a two-dimensional mesh in polar coordinates.

Formulation 2: Vortex-sheet formulation for twin jets
The vortex sheet model treats the shear layer as a boundary of infinitesimal width with a uniform streamwise velocity inside the jet and zero streamwise velocity outside (Lessen et al. 1965;Michalke 1970).The inviscid equations governing the linear instability waves, upon the introduction of the modal form (2.1) written in terms of the cylindrical coordindates, reduce to the Helmholtz equation for the disturbance pressure : This equation describes the disturbance pressure in the inner ( ) and outer ( ) regions to the vortex sheet, depending on the definition of : (2.6) In the twin-jet configuration (Morris 1990;Du 1993;Stavropoulos et al. 2021), the inner and outer solutions are written in forms consistent with the separation in even/odd families as ( 1 ) cos( 1 ) + ˆ ( 1 ) sin( 1), (2.7) where and are the modified Bessel functions of first and second kind.Each line of (2.8) corresponds to one of the four families.
The inner and outer solutions are matched at the ideally-expanded diameter , that is imposed to be equal to in this work.Matching conditions impose the continuity of the pressure and displacement across the vortex sheet: . (2.10) In order to impose the matching conditions, the outer solution (2.8) is re-written in terms of the coordinates of a single jet using Graf's addition theorem (Abramowitz & Stegun 1964): ( (2.12) Combining equations (2.7)-(2.12)and collecting terms corresponding to the same family yields the dispersion relation for a twin-jet system as where correspond to the coeffients , , or depending on the solution family, is the Kronecker delta and Here, is equal to 0.5 if = 0 and 1 otherwise.The mode families are imposed through the factors and in equations (2.13) and (2.15) (see table 1).For large the first factor in brackets in equation (2.15) tends to zero leaving only those from equation (2.14), which recovers the dispersion relation for a single jet (Lessen et al. 1965;Towne et al. 2017).As opposed to the case of a single jet, equation (2.13) couples all the azimuthal numbers.For the practical solution, the summatory is truncated at a finite .The results presented in this work are computed using = 5, but their convergence has been checked with a maximum of 10.

Finite-thickness single and twin-jet mean flows
The calculations of finite-thickness single and twin jets in this work employ analytic mean flows for the sake of reproducibility.The analytical velocity profile proposed by Michalke (1970Michalke ( , 1984) ) is used: where is the radial coordinate measured from the jet centre for a single jet and is the nozzle radius, = /2.A tailored mean flow field is constructed for twin jet configurations by combining the flow fields corresponding to two isolated round jets of the form (3.1) aligned with the −direction and with centres placed symmetrically at the coordinates ( , ) = (± /2, 0), as shown in figure 1.This is a good approximation of the twin-jet mean flow in the region immediate downstream of the nozzles, known as the converging region (Okamoto et al. 1985;Moustafa 1994), where the shear layers of the two jets are still well separated, and has been used in the past in the linear stability calculations by Morris (1990)  A zero mean pressure gradient is assumed and the Crocco-Busemann relation particularised for the temperature ratio , is used to determine the mean temperature ¯ .All the quantities in this expression are dimensionless, as explained in section 2. The mean density field ¯ is obtained from the state equation.
Following the parallel-flow assumption, cross-stream mean velocity components are neglected.
For the finite-thickness calculations herein, the jet Mach number is set at = 1.5 and = 1 and the parameter / = 12.5.This value of / is representative of the thin shear layer in the first diameter from the nozzle (Crighton & Gaster 1976;Michalke 1984).The Reynolds number = 5 × 10 4 is used in the calculations but identical results are recovered for = 10 5 , showing that the effect of viscosity is negligible.

The eigenspectra of single and twin jets
This section discusses the results of the LST analysis based on cartesian coordinates, described in section 2.2, when applied to single and twin jets.

Single round jets
When the local stability analysis is applied to single round jets, the axial symmetry of the mean flow allows for the introduction of azimuthal Fourier modes of the form exp(i ) (e.g.Gill 1965;Michalke 1984).This is usually exploited to reduce the eigenvalue problem to a onedimensional one, dependent on the radial coordinate alone.For each , the stability analysis recovers different families of eigensolutions (Tam & Hu 1989;Rodríguez et al. 2013), but only one discrete eigenmode is associated with the Kelvin-Helmholtz (K-H) instability.
The present approach is based on Cartesian coordinates, discretising only the first quadrant of the ( , )−plane and imposing symmetric/anti-symmetric conditions at the = 0 and = 0 planes.As a consequence, it does not isolate the individual azimuthal eigenmodes.Conversely, the imposed symmetries separate the eigenmodes in four families: SS, AS, SA and AA, as shown in table 1 and discussed elsewhere (Morris 1990;Rodríguez 2021;Rodríguez et al. 2021).When this approach is applied to a single jet centred at the origin of coordinates, each eigenspectrum contains a number of discrete K-H eigenmodes, each one corresponding to one .Figure 2  eigenvalues appear repeated in two families.In particular, the eigenmodes corresponding to = 1 are recovered in families AS and SA; in the following these eigenmodes are referred to as AS1 and SA1, respectively.Their pressure eigenfunctions are shown in figure 3( , ).The three-dimensional pressure fields are reconstructed following (2.1), with the eigenfunction ( , ) normalised.To aid in the visualisation, an arbitrary streamwise domain of length 10 is used and the pressure field amplitude is multiplied by exp(− ) to cancel the spatial growth, so that the amplitude remains constant.The pressure vanishes simultaneously at one of the coordinate axes (different for SA and AS) for the real and imaginary components.This behaviour is enforced here by the symmetry/anti-symmetry conditions for each family, but a computation using the domain = arctan(Im( )/Re( )), extracted at a cylinder of radius 0.75 aligned with the jet axis.The azimuthal angle 1 is measured from the positive axis as shown in figure 1.Only the subdomain 1 ∈ [0, ] is shown, but the omitted region can be inferred from the symmetries.Mode SA1 presents a constant phase for − /2 < 1 < /2 with is shifted an angle for other values of 1 .Mode AS1 shows a similar behaviour but with the phase shift at 1 = .Each of these eigenfunctions describes a jet flapping motion.However, because the eigenvalues corresponding to the modes AS1 and SA1 are identical ( 1 = 1 ), a linear combination of them is also a valid solution to the general problem.With appropriately chosen amplitude coefficients for each eigenmode, their linear combination recovers the features of an helical mode.This is illustrated in figure 3( , ), in which the pressure fields of modes AS1 and SA1 are added as ( , , , ) = ˆ SA1 ( , ) e i( SA1 − ) + i ˆ AS1 ( , ) e i( AS1 − ) + c.c.. (3. 3) The helical nature of the resulting mode is manisfest in the linear dependence of the phase on the azimuthal angle 1 .

Twin round jets
For the twin-jet configuration, the jet axes are not located at the origin of coordinates; while the nomenclature used for the single jet remains valid, the four solution families do not separate the odd and even modes.Instead, each family contains one eigenmode for each of the > 0 modes.The symmetry of the mean flow about the mid-plane leads to the recovery of the = 0 mode only in the SS and SA families.Representative eigenspectra for a twin jet configuration with separation / = 2.2 at = 0.3 is shown in figure 2( ), to be compared with the single jet eigenspectra in figure 2( ).
The interaction between the fluctuation fields of the two jets breaks the azimuthal symmetry of the eigenmodes and the same mode presents different eigenvalues , one for each symmetry family.This shift of the eigenvalues with respect to the single jet ones is stronger for = 0 and gradually reduces for increasing .This behaviour is explained by the asymptotic decay of the single-jet eigenfunctions as 1964).For a fixed frequency , the axial wavenumber increases with increasing , and hence leads to a faster radial decay of the fluctuations.Consequently, the amplitude of the pressure field associated with one jet that reaches the other jet is reduced.For the same reason, the magnitude of the eigenvalue shift is also inversely proportional to the jet separation / , as shown next.
In what follows attention is focused on the = 1 modes, which present the largest growth rates for most of the conditions considered.These modes are denoted by SS1, SA1, AS1 and AA1, where the two letters correspond to the symmetry family and the number to the corresponding azimuthal number.Figure 4 shows the eigenvalues corresponding to each of the four families as a function of the jet separation.While the magnitude of the eigenvalue shift grows continuously with decreasing / , the dependence of the real and imaginary parts of the eigenvalues is non-monotonic and different for each family.The shift is stronger for modes SS1 and SA1, which correspond respectivaly to in-phase (varicose) and counter-phase (sinuous) lateral flapping oscillations in the plane containing the jet axes.For jet separation above / ≈ 2, the real part of the wavenumber is reduced (equivalently, the phase velocity is increased) for SA1, while it is increased for SS1.This behaviour is reversed for very small jet separations.The flapping modes in the vertical (perpendicular) plane follow a similar trend: the in-phase mode AS1 reduces its phase velocity resulting from the presence of the other jet and the counter-phase AA1 increases it.As shown in figure 4( ) and more clearly in the zoomed-in version 4( ), the spatial growth rate presents a similar behaviour for the two pairs of modes: (i) The in-phase flapping modes SS1 (blue) and AS1 (red) exhibit an increase of the growth rate (− ) as / is reduced from / → ∞, due to jet interactions.For decreasing separations, the growth rate eventually reaches a maximum, then equals that of the single jet and then becomes smaller, i.e. close jet proximity stabilises the in-phase flapping modes.The jet separations for the destabilisation/stabilisation inversion are different for the two modes: / ≈ 2.8 for SS1 and 2.5 for AS1.The variation of the growth rate is of bigger amplitude for the varicose lateral flapping mode SS1.
(ii) The counter-phase modes SA1 (yellow) and AA1 (black) present opposite trends to the in-phase modes.These modes have the largest growth rates for small jet separations and, in particular, mode SA1 corresponding to sinuous lateral flapping is the most amplified one.
The eigenvalue shift resulting from the coupling of the unsteady twin-jet pressure fields has two inter-connected consequences.First, some of the = 1 eigenmodes become dominant over the others for a given combination of parameters ( / , , but also and , as shown later), suggesting that a preferred mode of flapping oscillation should emerge, as opposed to the helical motion of a single jet.From the results of figure 4, mode SS1 (varicose lateral flapping) should be expected for jet separation above / ≈ 3 and = 0.3.For / < 3, mode SA1 (sinuous lateral flapping) becomes gradually dominant.The second consequence of the eigenvalue shift is that = 1 helical modes cannot be recovered as the linear combination of different = 1 modes.Owing to the exponential dependence on (see equation 3.3), even small differences in the spatial growth rate result in a fast dominance of the most unstable flapping eigenmode.This is illustrated in figure 5 for a jet separation of / = 2.2.As was done for figure 3, the streamwise amplification of the eigenfunctions is re-scaled with the growth rate of the most unstable one (SA1 in this case).As a consequence, the least amplified mode AA1 has its relative amplitude reduced along the streamwise direction.A phase mismatch also appears due to the small differences in .The combination of the two modes is shown in figure 5( , ): the pressure field behaves initially as two anti-symmetric counter-rotating helical modes, but soon evolves towards the sinuous lateral flapping of mode SA1.

Preferred mode of oscillation for finite-thickness twin jets
The results in the previous section suggest a strong dependence of the eigenvalue shift for all modes with frequency and jet separation.Figure 6 shows the results of a parametric study of the = 1 modes for varying between 0.1 and 0.6 and / varying between 1.8 and 5, extending the results of figure 4 for the twin jets with = 1.5, = 1, and / = 12.5.Figure 6( ) shows the parametric regions in which each mode dominates, i.e. its spatial growth rate is the largest among the = 1 modes.Incidentally, it is noted that the growth rate of all = 1 is always larger than that of the = 0 modes, as expected for high jets with relatively thin shear layer (Michalke 1984).The alternating dependence of the leading mode with / described before is also obtained with , and the region of preferrence of the modes forms "stripes" whose boundaries take approximately constant values of the product × / .
Figure 6( ) provides no information on the relative growth rates of the leading mode with respect the others, or with respect to eigenmodes corresponding to a single jet. Figure 6( ) shows, for the leading eigenmode at each ( , / ), the relative change in the growth rate resulting from the jet interaction, quantified as Note that the dispersion relation for twin jets ( , / ) recovers that of a single jet in a smooth manner as / → ∞.Each colour level in figure 6( ) corresponds to a relative increase of 1% of Δ .The maximum colour level shown is 10% but Δ reaches larger values.As anticipated, the destabilisation of the leading eigenmode is inversely proportional to jet separation.For a fixed jet separation, the envelope of |Δ | over all modes is also inversely proportional to , similarly to the dependence on / shown in figure 6(d).Consequently, as the product × / increases, Δ approaches zero.This happens in this case for regions dominated by the modes SS1 and AS1. Figure 6( ) shows again Δ , but colour coded to highlight the preferred oscillation mode.For the jet conditions analysed herein ( = 1.5, = 1, / = 12.5), these results predict the dominance of mode SA1 leading to lateral sinuous oscillations for most jet separations and Strouhal numbers.The white regions in figure 6( ) correspond to conditions in which the jet interactions are weak, and hence do not lead to a substantial eigenvalue shift; the oscillation modes converge towards their behaviour for single jets for which there is no preferred mode.In this case, the jet may follow either helical of flapping dynamics, which may be selected by other characteristics of the turbulent flow.
The results presented so far consider finite-thickness jets and are obtained using the linear stability analysis formulation described in section 2.1.While those calculations are performed using a personal computer, the parametric study leading to figure 6 requires the numerical solution of over 10 000 individual matrix eigenvalue problems, which becomes impractical for more comprehensive studies including the effect of the jet Mach number and the temperature ratio.In the following section, such a study is performed using the (zero-thickness) twin-jet shear layer model presented in section 2.2.

Dependence of the preferred mode of oscillation on Mach number and temperature ratio
This section presents a parametric study of the effect of the jet Mach number and temperature ratio on the preferred oscillation mode, using the vortex-sheet model for twin round jets.The description of the eigenmode families, the eigenvalue shift and their qualitative dependence on jet separation and Strouhal number given for finite-thickness jets holds for the vortex-sheet model.For each ( , ), the evolution of the four = 1 eigenmodes with jet separation and Strouhal number is monitored in a manner analogous to that described in the previous section.The reduced computational cost of the vortex-sheet model computations allows for a greater region of interest and to significantly refine the discretisation of the ( , / ) space.The analysis for one ( , ) case involves the calculation of nearly 1.5 million eigenvalues, but takes less than an hour on a modern personal computer.
Figure 7 shows the leading oscillation mode for twin jets with = 1 and between 1 and 2. The overall picture reproduces that already described for the finite-thickness = 1.5 case: (i) the regions in which each particular mode dominates in the ( , / ) space form alternating stripes; (ii) the maximum eigenvalue shift, computed following (3.4), occurs for the smallest jet separation within each strip; (iii) as the jet separation is reduced, the Strouhal number for the maximum shift of the growth rate is increased.For = 1 (figure 7( )) and / 3, mode SS1 dominates for 0.2.For larger , modes SA1 and AA1 become dominant but their eigenvalue shift is comparatively weak.Increasing the jet Mach number at constant temperature ratio has two effects.First, the boundaries of the mode stripes are gradually displaced towards lower and / values.Second, the shift of the growth rates is increased.The latter implies that the interaction between the fluctuation fields of the two jets intensifies with .The combination of the two effects leads to the gradual appearance of a new stripe of mode SS1 at higher while the corresponding stripe at the lower region reduces.Simultaneously, the stripe of mode SA1 is displaced to lower Strouhal numbers while its growth rate shift is increased.For above 1.5, mode SA1 (sinous lateral oscillations) has the largest growth rate shift.
For above 1.6 (figures 7( − )), the results for mode SA1 present an anomalous behaviour in the limit of lower separation and Strouhal number, that is clearly visible in the figures.Instead of following the general trend, boundary of the stripe forms a horizontal segment (constant ) for / below a threshold value increases with .A preliminary study revealed that this anomaly is due to the appearance of a saddle point in the evolution of mode SA1, which pinches with mode SA2 at a certain ( , / ) combination.This saddle point is not associated with an absolute instability, as it involves two downstream-propagating K-H eigenmodes.A deep study regarding this issue is currently underway but is beyond the scope of this paper.
Figure 8 shows the parametric study for twin jets with fixed jet Mach number = 1.5 and temperature ratio increasing from = 1 to 4. The effect of increasing the jet temperature follows the same trends as increasing the jet Mach number: the mode stripes move towards lower and / , reducing the region dominated by mode SS1 (dominant for = 1) and increasing the relative importance of mode SA1 with increasing temperature.Overall, the impact of the temperature ratio on the eigenmode shift is weaker than that of the jet Mach number, but still mode SA1 attains growth rates nearly 10% larger that those for the equivalent single jet.

Conclusions
Linear stability theory is often used to model the large-scale flow structures in the turbulent mixing region and near pressure field of high-speed jets.For perfectly-expanded single round jets, these models predict the dominance of = 1 helical modes for relatively low frequencies and thin shear layers, in agreement with empirical data.
For twin-jet configurations, the interaction between their flutuation fields results in coupling.The symmetries present in the geometry lead to a separation of the modes into four families, corresponding to the even/odd behaviour of the pressure field about the symmetry plane and the plane containing the jet axes.The division into families is inherent to the twin-jet geometry and not just an artifact of the problem formulation: while it is exploited here to simplify the analyses, a formulation in which the symmetries/anti-symmetries are not imposed recovers the same families (Rodríguez et al. 2018).The fluctuation field associated with = 1 modes for twin jets does not correspond to helical oscillations, but to flapping oscillations of the jet.
The coupling between the jets is quantified by monitoring the shift of the complex wavenumber (eigenvalue of the spatial linear stability problem) with respect to that of the single jet, Δ .The absolute value of Δ is different for each oscillation mode, but for all of them is inversely proportional to both the jet separation / and the Strouhal number .However, the real and imaginary parts of Δ do not present a monotonous dependence on or / .This leads to regions in the ( , / ) space in which different modes of oscillation dominate over others.A parametric study is presented in which the jet Mach number is varied between 1 and 2, and the jet temperature ratio is varied between 1 and 4. Increasing independently both and is found to augment the jet coupling and modify the ( , / ) map of the preferred oscillation mode.Present results predict that only two oscillation modes are likely to be observed when jet coupling is relevant, namely in-phase (varicose) or counter-phase (sinuous) flapping oscillations on the plane containing the nozzles (modes SS1 and SA1,respectively) It should be noted that while flapping modes are the natural mode of oscillation of twin jets, a superposition of them may still generate some helical motions in regions close to the nozzle.However, owing to the differences between the growth rates of each mode and their exponential amplification in the first few diameters, it can be expected that flapping motions will dominate the flow as it develops downstream.Hence, the eventual helical oscillations in twin-jet systems would be a spatially transient feature of the flow.
The results presented in this work consider perfectly expanded jets.The preference of certain oscillation modes over others is predicted based on the larger growth rate of the corresponding Kelvin-Helmholtz instability eigenmode.It may be expected that these results would be applicable also to shock-containing jets, at least to some extent.The preference of shock-containing twin jets to present flapping motions at conditions in which a single jet presents helical motions has been shown by different authors (Seiner et al. 1988;Wlezien 1989;Kuo et al. 2017;Knast et al. 2018;Bell et al. 2018Bell et al. , 2021)).However, one should be cautious when extending the present conclusions to shock-containing jets.As these jets are subjected to the resonance loop associated with screech (Powell 1953;Edgington-Mitchell 2019;Nogueira et al. 2022b) (2022), the closure mechanism of screech is strongly dependent on the wavenumber of the waves involved in resonance, and less dependent on the growth rate of the Kelvin-Helmholtz instability.However, as the wavenumber of the Kelvin-Helmholtz eigenmode is one of the driving components of this phenomenon, the present results remain relevant: they may be used as input for screech models in twin jets and help explain the preferred flapping mode in shock-containing jets as well.
Twin-jet configuration and geometry, showing the different coordinate systems employed.
( )   shows the eigenspectra for a single jet at = 0.3.The different symbols correspond to one of the solution families.Except for = 0, which is only recovered in family SS, all other

F 3 .
Pressure eigenfunctions corresponding to a single jet at = 1.5, = 1, / = 12.5 and = 0.3.The streamwise dependence is obtained from the corresponding eigenvalue, eliminating the spatial growth.Left column: iso-contours of the real pressure component.The grey circle shows the nozzle circumference.Right column: phase angle at a cylinder of radius 0.75 .( , ) SA1; ( , ) AS1; ( , ) Linear combination of SA1 and AS1 to produce the = 1 helical mode.

F 5 .
Pressure eigenfunctions corresponding to twin jets with / = 2.2 at = 1.5, = 1, / = 12.5 and = 0.3.The streamwise dependence is obtained from the corresponding eigenvalue, eliminating the spatial growth corresponding to SA1.Left column: iso-contours of the real pressure component.The grey circles shows the nozzle circumference.Right column: phase angle at a cylinder of radius 0.75 centred on one jet.( ) SA1; ( ) AA1; ( ) Linear combination of SA1 and AA1.
Preferred oscillation mode for twin jets as a function of the jet separation and Strouhal number.= 1.5, = 1, / = 12.5.( ) Leading eigenmode.( ) Relative increase of the growth rate with respect to the single jet.( ) Same as ( ), but colour-coded to show the leading eigenmode.The same colour-coding of figure 4 is used: Blue: SS1; Red: AS1; Yellow: SA1; Black: AA1.
. The azimuthal dependence for each family is also shown.The last two columns show the values of and appearing in the vortex sheet model.