1. Introduction
Turbulent flows are ubiquitous in nature and industrial applications, although their interpretation is challenging due to their inherent complexity. A common strategy to tackle their study relies on coherent structures – organised regions in space with similar features that persist in time – which provide a low-dimensional representation of otherwise high-dimensional turbulent flows. This approach has been exploited for more than half a century, with some early works such as those by Kline et al. (Reference Kline, Reynolds, Schraubt and Runstadlers1967) who revealed streaky patterns in turbulent boundary layers, or Cros & Champagne (Reference Cros and Champagne1971) and Brown & Roshko (Reference Brown and Roshko1974) who identified organised structures in turbulent jets and mixing layers, respectively. More recently, the dynamics of coherent structures has been linked to relevant engineering problems, such as drag reduction in wall-bounded turbulent flows (Jiménez Reference Jiménez2018) and noise generation in high-speed turbulent jets (Jordan & Colonius Reference Jordan and Colonius2013).
Here, we focus on large-scale coherent structures, a few of which are sufficient to picture the dominant flow dynamics. There are different ways of computing these large-scale coherent structures, with mathematical tools such as modal decompositions being a common data-driven approach for their identification (Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017). These algorithms aim to describe the flow as a superposition of spatial modes, each accompanied by characteristic values that represent either their energy content or dynamic traits. Modal decompositions generally require large amounts of data, and can be used to compress high-dimensional data (hundreds thousands of grid points in a numerical simulation) into more interpretable low-dimensional forms. This simpler description of turbulent flows is also suitable for reduced-order models, where modal decompositions can be combined with Galerkin methods (Rowley & Dawson Reference Rowley and Dawson2017), run in parallel with numerical simulations (Amor et al. Reference Amor, Schlatter, Vinuesa and Le Clainche2023) or integrated into deep learning models (Abadía-Heredia et al. Reference Abadía-Heredia, López-Martín, Carro, Arribas, Pérez and Le Clainche2022; Le Clainche, Rosti & Brandt Reference Le Clainche, Rosti and Brandt2022) to accelerate computations.
In this work, we employ modal decompositions to identify coherent structures in turbulent planar jets. In particular, we use higher-order dynamic mode decomposition (HODMD) (Le Clainche & Vega Reference Le Clainche and Vega2017). Similarly to standard dynamic mode decomposition (DMD) (Schmid Reference Schmid2010), HODMD extracts the dominant nonlinear dynamics directly from the data. The HODMD (and DMD) represents the flow as a linear superposition of so-called modes, each corresponding to coherent structures that oscillate either in time, space or in time and space. The key innovation of HODMD lies in its use of delay embeddings, which augment the original data by incorporating time-delayed snapshots (Takens Reference Takens1981). This improves the spectral resolution of the data and filters small-amplitude frequencies (Le Clainche et al. Reference Le Clainche and Vega2017), thus making the method more suitable for analysing highly complex, nonlinear flows (Le Clainche et al. Reference Le Clainche, Izbassarov, Rosti, Brandt and Tammisola2020, Reference Le Clainche, Rosti and Brandt2022). Here, we extend the application of HODMD to compute the coherent structures in Newtonian and viscoelastic turbulent planar jets. Specifically, we validate the method for high-Reynolds-number Newtonian planar jets, and we extend the work by Amor et al. (Reference Amor, Corrochano, Foggi Rota, Rosti and Le Clainche2024) for low-Reynolds-number viscoelastic planar jets.
Coherent structures in Newtonian turbulent jets are characterised by three distinct mechanisms. First, the Kelvin–Helmholtz instability causes wave packets – coherent axisymmetric structures – that grow linearly near the inlet (Crighton & Gaster Reference Crighton and Gaster1976; Cohen & Wygnanski Reference Cohen and Wygnanski1987; Suzuki & Colonius Reference Suzuki and Colonius2006; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013). A second mechanism, the Orr mechanism, occurs downstream of the potential core, where the flow response is rather nonlinear (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013; Tissot et al. Reference Tissot, Lajús, Cavalieri and Jordan2017; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018). In this region, Orr wave packets have greater gain than those associated with Kelvin–Helmholtz instability, specifically at low frequency and zero wavenumber (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018). The third response mechanism is the lift-up mechanism (Boronin, Healey & Sazhin Reference Boronin, Healey and Sazhin2013; Jiménez-González & Brancher Reference Jiménez-González and Brancher2017; Nogueira et al. Reference Nogueira, Cavalieri, Jordan and Jaunet2019). The lift-up mechanism and the associated coherent structures, often called streaks, are an important mechanism in wall-bounded (Brandt Reference Brandt2014) and planar free-shear flows (Pierrehumbertt & Widnall Reference Pierrehumbertt and Widnall1982; Bernal & Roshko Reference Bernal and Roshko1986; Metcalfe et al. Reference Metcalfe, Orszag, Brachet and Riley1987). In turbulent jets, streamwise vortices have a fundamental effect on the dynamics and statistical properties of the flow (Liepmann & Gharib Reference Liepmann and Gharib1992), since these streamwise vortices (and streaks) are the most energetic for non-zero wavenumbers (Citriniti & George Reference Citriniti and George2000), scaling inversely with the distance from the inlet (Jung, Gamard & George Reference Jung, Gamard and George2004; Tinney, Glauser & Ukeiley Reference Tinney, Glauser and Ukeiley2008).
A comparable characterisation of coherent structures in viscoelastic jets is still lacking. Viscoelasticity, for instance introduced by adding flexible polymers to a Newtonian solvent, significantly changes the flow dynamics and, as a consequence, it is expected to also modify the coherent structures. In viscoelastic turbulent flows, polymers interact with the energy cascade by extracting energy from the large-scale eddies, that they either return to the flow at smaller scales or dissipate through polymer diffusion (Valente, Da Silva & Pinho Reference Valente, Da Silva and Pinho2014). Even though their influence is more evident at small scales, their range of influence, identified by changes in the scaling exponent of the energy spectrum, depends on the polymer elasticity (Rosti, Perlekar & Mitra Reference Rosti, Perlekar and Mitra2023) and inertia (Singh & Rosti Reference Singh and Rosti2025). Polymers can also destabilise laminar flows through purely elastic instabilities (Larson Reference Larson1992; Shaqfeh Reference Shaqfeh1996), leading to elastic turbulence (Groisman & Steinberg Reference Groisman and Steinberg2000; Singh et al. Reference Singh, Perlekar, Mitra and Rosti2024). These have been extensively studied in curvilinear flow configurations, such as the flow between concentric cylinders (Larson, Shaqfeh & Muller Reference Larson, Shaqfeh and Muller1990; Groisman & Steinberg Reference Groisman and Steinberg1998), counter-rotating parallel disks (McKinley, Byars & Armstrong Reference McKinley, Byars and Armstrong1991; Öztekin & Brown Reference Öztekin and Brown1993) and in curved channels (Groisman & Steinberg Reference Groisman and Steinberg2001a , Reference Groisman and Steinberg2004). However, curvature is not a prerequisite for their occurrence, and elastic instabilities can occur in straight shear flows too (Pan et al. Reference Pan, Morozov, Wagner and Arratia2013; Qin et al. Reference Qin, Salipante, Hudson and Arratia2017; Jha & Steinberg Reference Jha and Steinberg2021; Lellep, Linkmann & Morozov Reference Lellep, Linkmann and Morozov2024; Rota et al. Reference Rota, Amor, Le Clainche and Rosti2024). Recently, elastic turbulence has been found also hiding in the smallest scales of polymeric turbulence at large Reynolds number (Garg & Rosti Reference Garg and Rosti2025).
In jet flows, elastic instabilities appear at markedly lower Reynolds number compared with classic Newtonian jets. Yamani et al. (Reference Yamani, Keshavarz, Raj, Zaki, McKinley and Bischofberger2021) observed in round jets the transition to elasto-inertial turbulence – when inertial and elastic effects are of the same order of magnitude (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013). The elastic turbulent state can be sustained at much lower Reynolds number if elasticity is sufficiently large (Soligo & Rosti Reference Soligo and Rosti2023). In the latter case, the destabilising role of elasticity remains unclear. While elasticity stabilises the sinuous (antisymmetric) mode and partly stabilises the varicose (symmetric) mode at high Reynolds number (Rallison & Hinch Reference Rallison and Hinch1995), elasticity is rather destabilising at moderate Reynolds numbers, although this effect has a non-monotonic trend: increasing elasticity first destabilises the flow, although it stabilises if elasticity is sufficiently large (Ray & Zaki Reference Ray and Zaki2015), as also confirmed experimentally (Yamani et al. Reference Yamani, Keshavarz, Raj, Zaki, McKinley and Bischofberger2021, Reference Yamani, Raj, Zaki, McKinley and Bischofberger2023). This effect may be induced by the competing influence of elasticity in separated regimes: in the linear regime, elasticity enhances the instability of the jet (Ray & Zaki Reference Ray and Zaki2015), but it inhibits the roll-up process in the nonlinear regime, which yields overall flow stabilisation (Kumar & Homsy Reference Kumar and Homsy1999; Guimarães, Pinho & da Silva Reference Guimarães, Pinho and da Silva2023). In the low-Reynolds-number limit of elastic turbulence, polymers alone sustain the turbulence, when inertial effects are negligible (Groisman & Steinberg Reference Groisman and Steinberg2000).
To address this gap and clarify the mechanism, we apply nonlinear DMD (HODMD) (i) to compare the global coherent structures in Newtonian and viscoelastic turbulent planar jets and (ii) to examine the influence of polymer elasticity in the near-field region. Data-driven analysis is well suited for investigating global instabilities in highly nonlinear systems, where traditional stability analysis is often intractable. For instance, the spatial decomposition based on HODMD (or DMD) can represent the flow as growing and/or decaying time-dependent perturbations that propagate upstream or downstream from their point of origin. This approach was successfully applied by Yamani et al. (Reference Yamani, Raj, Zaki, McKinley and Bischofberger2023), who employed local DMD analysis to characterise elasto-inertial instabilities in viscoelastic planar jets. In our study, we decompose the nonlinear dynamics through a spatio-temporal analysis, in which the sequential application of HODMD yields a representation of the flow as a superposition of travelling waves. This methodology, namely spatio-temporal Koopman decomposition (Le Clainche & Vega Reference Le Clainche and Vega2018), has previously proven effective in capturing the unstable modes of the flow past a circular cylinder at low Reynolds number (Le Clainche et al. Reference Le Clainche and Vega2018), as well as the modes associated with streak breakdown in the turbulent channel flow of an elastoviscoplastic fluid (Le Clainche et al. Reference Le Clainche, Izbassarov, Rosti, Brandt and Tammisola2020).
The article is organised as follows. Section 2 introduces the direct numerical simulations and the data-driven analysis employed in this work. Section 3 compares the coherent structures in Newtonian and viscoelastic turbulent planar jets, with a detailed analysis of the near-field structures in the viscoelastic jet in § 4. Section 5 focuses on the polymer field in the viscoelastic jet, and finally § 6 summarises the main conclusions.
2. Methods
2.1. Direct numerical simulations
Data are obtained by means of direct numerical simulations. We employ our in-house solver Fujin (Rosti Reference Rosti2026) (https://www.oist.jp/research/research-units/cffu/fujin) for simulating the Newtonian and viscoelastic turbulent planar jets. Most of the data considered here are taken from previously published results (Soligo & Rosti Reference Soligo and Rosti2023; Soligo, Chiarini & Rosti Reference Soligo, Chiarini and Rosti2025).
The problem is governed by the incompressible Navier–Stokes equations
\begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0, \end{equation}
\begin{equation} \rho \left (\partial _t \boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u} \boldsymbol{u} \right ) = - \boldsymbol{\nabla }\!p + \mu _s {\nabla} ^2 \boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{\cdot }\tau , \end{equation}
where
$\boldsymbol{u}$
and
$p$
are the velocity and pressure fields,
$\rho$
is the density,
$\mu _s$
is the dynamic viscosity of the solvent and
$\tau$
is the non-Newtonian stress tensor. In the Newtonian jet,
$\tau = 0$
and
$\mu _s$
is equivalent to the viscosity of the fluid (
$\mu _0 = \mu _s$
). On the other hand, the viscoelastic jet requires an extra equation to model
$\tau$
. Here, we use the Oldroyd-B model
\begin{equation} \lambda \overset {\boldsymbol{\nabla }}{\tau } + \tau = \mu _{\!p} \big( \boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^{\top } \big), \end{equation}
where
$\lambda$
is the relaxation time of the polymer, i.e. the time required by the polymer to reach equilibrium if perturbed by an external forcing, and
$\mu _{\!p}$
is the dynamic viscosity of the polymer, with
$\mu _0 = \mu _s + \mu _{\!p}$
being the total viscosity. Here,
$\overset {\boldsymbol{\nabla }}{\tau }$
is the upper-convective derivative, defined as
\begin{equation} \overset {\boldsymbol{\nabla }}{\tau } = \partial _t \tau + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\tau - \big( \boldsymbol{\nabla }\boldsymbol{u}^{\top } \boldsymbol{\cdot }\tau + \tau \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} \big ). \end{equation}
Note that the model describes a purely viscoelastic fluid without shear-thinning viscosity. The non-Newtonian stress
$\tau$
can be written in terms of the conformation tensor
$\boldsymbol{C}$
, a second-order, positive-definite tensor that represents the average value of the end-to-end distance of the polymers, i.e.
$\tau ( \boldsymbol{C} ) = \mu _{\!p} ( \boldsymbol{C} - \boldsymbol{I} )/ \lambda$
, with
$\boldsymbol{I}$
the tensorial identity.
The equations are discretised on a staggered, uniform, Cartesian grid, with
$x$
being the streamwise,
$y$
the jet-normal and
$z$
the spanwise direction. Scalar variables, namely pressure, density, viscosity and extra stresses, are stored at the centre of the cells, and velocity at the faces. The spatial derivatives are approximated using second-order, central finite differences, while time integration is performed using a second-order explicit Adams–Bashforth scheme, coupled with a fractional step method (Kim & Moin Reference Kim and Moin1985) to solve the pressure coupling. We adopt a matrix-logarithm formulation (Fattal & Kupferman Reference Fattal and Kupferman2004; Hulsen, Fattal & Kupferman Reference Hulsen, Fattal and Kupferman2005) to overcome the numerical instability at high polymer elasticity (Keunings Reference Keunings1986), and we use a fifth-order weighted essentially non-oscillatory scheme (Shu Reference Shu2009; Sugiyama et al. Reference Sugiyama, Ii, Takeuchi, Takagi and Matsumoto2011) to solve the advection term in the upper-convective derivative.
The fluid is injected through a plane slit with height
$2h$
and length
$L_z$
from the left side of the domain box, that is filled with the same fluid as that injected through the plane slit. We impose on the left boundary (
$x = 0$
) the no-slip and no-penetration boundary conditions, except for the inlet slit (centred at
$y = 0$
), where a top-hat profile is imposed to the streamwise velocity (Da Silva & Métais Reference Da Silva and Métais2002; Stanley, Sarkar & Mellado Reference Stanley, Sarkar and Mellado2002)
\begin{equation} \frac {u}{U} = \frac {1 + \alpha }{2} + \frac {1 - \alpha }{2} \textrm {tanh} \left [ \frac {h}{4 \vartheta } \left ( \frac {2 |y|}{h} - 1 \right ) \right ]\!, \end{equation}
where
$\vartheta$
is the shear layer momentum thickness, and
$\alpha = 0.1$
in the Newtonian case and
$\alpha = 1$
in the viscoelastic case; a small co-flow of
$0.1 U$
is added in the Newtonian jet. This choice is made to minimise the effect of the upper and lower boundaries in the Newtonian jet, where, as we describe below, the vertical length of the computational box is significantly smaller compared with that for the viscoelastic jet. Simulating the Newtonian jet in a similar domain to the viscoelastic jet would be computationally unfeasible, thus the different choice of co-flow. Free-slip and no-penetration boundary conditions are imposed in the lower and upper boundaries (
$y = \pm L_y / 2$
), periodic boundary conditions are used in the spanwise direction (
$z = 0$
,
$z = L_z$
) and a non-reflective outflow boundary condition (Orlanski Reference Orlanski1976) is enforced at the outlet (
$x = L_x$
).
The Reynolds number is based on the half-slit height
$h$
and the inlet velocity
$U$
, and is defined as
$\textit{Re} = \textit{Uh} / \nu _0$
, where
$\nu _0$
is the total kinematic viscosity of the fluid, equal to
$\mu _s/\rho$
in the Newtonian case and to
$ ( \mu _s +\mu _{\!p} )/\rho$
in the viscoelastic one. We consider a Newtonian jet at a high Reynolds number,
$\textit{Re} = 6750$
, and a viscoelastic jet at a low Reynolds number,
$\textit{Re} = 20$
. We choose a large enough Reynolds number to achieve fully developed turbulence in the Newtonian jet (Soligo et al. Reference Soligo, Chiarini and Rosti2025); contrarily, Newtonian jets at
$\textit{Re} = 20$
are laminar, so any observed unstable behaviour in the viscoelastic jet is attributed to the non-Newtonian property of the fluid (Soligo & Rosti Reference Soligo and Rosti2023). The viscoelastic jet requires two additional non-dimensional parameters. The first one is the Deborah number,
$De = \lambda U / h$
, that is set to
$100$
, such that the elasticity number, measured as
$E = \textit{De} / \textit{Re}$
, is greater than unity; in other words, turbulence in this regime is sustained solely by elasticity, with inertial effects being negligible. We denote this regime elastic turbulence even though
$\textit{Re}$
is not smaller than unity, since
$\textit{Re} \lesssim De$
and since the flow would be laminar in the absence of polymers (Groisman & Steinberg Reference Groisman and Steinberg2001b
). At the same
$\textit{Re}$
but one order of magnitude smaller
$De$
(
$De = 10$
), we have
$E \lesssim 1$
, and the resulting viscoelastic jet has properties similar to those of the elastic turbulent planar jet studied here, but with the flow relaminarising downstream after a certain distance from the inlet (Soligo & Rosti Reference Soligo and Rosti2023). Once fully developed, elastic turbulence shows similar qualitative properties notwithstanding the value of
$De$
(Singh et al. Reference Singh, Perlekar, Mitra and Rosti2024). Finally, the last parameter is the viscosity ratio
$\beta = \mu _s / \mu _0 = \mu _s / ( \mu _s + \mu _{\!p} )$
, that is set to
$0.98$
, thus indicating a dilute polymer solution.
The computational domain extends for
$140h \times 140h \times 35h$
in the streamwise, jet-normal and spanwise directions for the Newtonian case, and
$320h \times 480h \times 26.7h$
in the viscoelastic case, and each domain is discretised using
$5760 \times 5760 \times 1440$
and
$1440 \times 2340 \times 128$
grid points, respectively. The size of the two domains are different to ensure a minimal effect of the bounding box on the flow behaviour. Indeed, due to the difference of the Reynolds numbers, the computational box in the low-Reynolds-number viscoelastic jet must be longer and taller to contain a fully developed region of the flow, which develops further away from the inlet, compared with the Newtonian turbulent jet. The level of turbulence in the viscoelastic jet remains, however, much less intense than in the Newtonian jet: the Taylor Reynolds number in the Newtonian jet reaches
$\textit{Re}_{\lambda } \approx 120$
in the developed region beyond
$x/h \gtrsim 80$
, while for the viscoelastic one it is limited to
$\approx 4$
at the same distance. The chosen set-up is the results of various tests done with different box sizes, performed by Soligo & Rosti (Reference Soligo and Rosti2023) and Soligo et al. (Reference Soligo, Chiarini and Rosti2025), respectively.
2.2. Higher-order dynamic mode decomposition
Higher-order dynamic mode decomposition (Le Clainche & Vega Reference Le Clainche and Vega2017) is an extension of the standard DMD (Schmid Reference Schmid2010). Similarly to DMD, HODMD decomposes spatio-temporal data
$\boldsymbol{v} ( \boldsymbol{x}, t )$
into a finite set of
$M$
DMD modes, each associated with a frequency
$\omega _m$
and a growth rate
$\delta _m$
. In addition, HODMD computes an amplitude coefficient
$a_m$
that weights each spatial mode
$\boldsymbol{u}_m(\boldsymbol{x})$
in the reconstruction, quantifying the contribution of each mode to the overall flow dynamics. More precisely, HODMD (and DMD) reconstructs the data as follows:
\begin{equation} \boldsymbol{v} \left ( \boldsymbol{x}, t_k \right ) \simeq \sum _{m = 1}^{M} a_m \boldsymbol{u}_m \left ( \boldsymbol{x} \right ) e^{\left ( \omega _m i + \delta _m \right ) t_k}, \quad k = 1, \ldots , N_t. \end{equation}
To estimate (2.6), the data are arranged in matrix form
$\boldsymbol{V}_1^{N_t}$
, with each column containing a snapshot
$\boldsymbol{v}_k$
– an instantaneous field from an experiment or numerical simulation. This matrix has dimension
$I \times J \times N_t$
, where
$I$
is the number of flow features,
$J = N_x \times N_y \times N_z$
the number of spatial grid points and
$N_t$
the total number of snapshots – they should be spaced equally in time.
There are two fundamental differences between HODMD and the standard DMD method. First, HODMD introduces a pre-processing dimensionality reduction step. The snapshot matrix
$\boldsymbol{V}_1^{N_t}$
is projected onto a low-dimensional basis, yielding a compact representation that exploits redundancies and mitigates noise. Second, HODMD augments the projected data using delay embeddings (Takens Reference Takens1981). The combination of DMD with time-delayed observables, e.g. in the form of a Hankel matrix, has proven effective for approximating the spectral properties of the Koopman operator in nonlinear dynamical systems (Arbabi & Mezic Reference Arbabi and Mezic2017; Brunton et al. Reference Brunton, Brunton, Proctor, Kaiser and Kutz2017; Kamb et al. Reference Kamb, Kaiser, Brunton and Kutz2020). This implementation makes HODMD suitable for extracting spectral information from temporally broadband data, which enables carrying out of the standard DMD method in highly complex nonlinear flows.
In what follows, we outline the steps comprising HODMD. For further details on DMD, the reader is referred to Schmid (Reference Schmid2010).
2.2.1. Step 1: dimensionality reduction via singular value decomposition
Truncated singular value decomposition (SVD) is performed to the snapshot matrix
\begin{equation} \boldsymbol{V}_1^{N_t} \simeq \boldsymbol{U}\! \boldsymbol{\varSigma } \boldsymbol{T}^{\top }, \end{equation}
where the diagonal matrix
$\boldsymbol{\varSigma }$
contains the singular values
$\sigma _1 \gt \sigma _2 \gt \ldots \gt \sigma _N$
, and
$\boldsymbol{U}^{\top } \boldsymbol{U} = \boldsymbol{T}^{\top } \boldsymbol{T}$
are
$N \times N$
unit matrices. Based on a threshold
$\varepsilon _1$
, the spatial dimension
$I \times J$
of the original data is reduced to a set of linearly independent vectors of dimension
$N$
, with
$N \ll I \times J$
the spatial complexity of the reduced data. The reduced snapshot matrix is defined as
\begin{equation} \hat {\boldsymbol{V}}_1^{N_t} = \boldsymbol{\varSigma } \boldsymbol{T}^{\top }, \end{equation}
and has dimension
$N \times N_t$
. The reduced snapshot matrix is computationally more tractable than the original one, thus making possible the implementation of the delay embeddings, in particular if data are very high-dimensional, which is the usual case in large-scale problems.
In this work, we implement high-order SVD (HOSVD) (Tucker Reference Tucker1966) instead of standard SVD. The HOSVD is better suited for tensorial data, as it captures the relationships among different dimensions. Truncation can be performed separately for each dimension, allowing the complexity to be adjusted independently. This is advantageous for our case, where planar jets exhibit greater complexity in the streamwise and jet-normal directions than in the spanwise one. However, this flexibility comes at the cost of increased parametrisation: the number of thresholds increases from one (uniform truncation in space and time) to four (three in space and one in time). To simplify this, we adopt either a single threshold
$\varepsilon _{1}$
applied to all dimensions, or two thresholds,
$\varepsilon _{1s}$
in space, and
$\varepsilon _{1t}$
in time. Although HOSVD increases the computational cost, its use provides a more optimal implementation of HODMD for multi-dimensional data (Le Clainche et al. Reference Le Clainche and Vega2017).
2.2.2. Step 2: time-delay embedding
Next, HODMD introduces the following assumption:
\begin{equation} \hat {\boldsymbol{V}}_{d+1}^{N_t} \simeq \hat {\boldsymbol{R}}_1 \hat {\boldsymbol{V}}_1^{N_t - d} + \hat {\boldsymbol{R}}_2 \hat {\boldsymbol{V}}_2^{N_t - d + 1} + \cdots + \hat {\boldsymbol{R}}_d \hat {\boldsymbol{V}}_d^{N_t - 1}, \end{equation}
where the operator
$\hat {\boldsymbol{R}}$
maps the temporal evolution of the reduced observables onto an infinite-dimensional linear space (Rowley et al. Reference Rowley, Mezic, Bagheri, Schlatter and Henningson2009; Mezic Reference Mezic2013). Standard DMD exploits this concept and HODMD refines it by incorporating time-lagged snapshots that generalise the original Koopman assumption from DMD to the high-order expression in (2.9). Figure 1 illustrates this process, where a window of length
$N_t - d$
slides through
$\hat {\boldsymbol{V}}$
. Note that HODMD reduces to standard DMD when
$d = 1$
(see also in (2.9)).
Sketch of the delay embedding. A window of length
$N_t - d$
(blue outline) moves through the reduced snapshot matrix, producing the delayed matrix
$\hat {\boldsymbol{V}}_{d+1}^{N_t}$
(red box), which is expressed as a linear superposition of the
$d$
preceding reduced snapshot matrices.
Equation (2.9) can be recast as
\begin{equation} \tilde {\boldsymbol{V}}_2^{N_t - d + 1} = \tilde {\boldsymbol{R}} \tilde {\boldsymbol{V}}_1^{N_t - d}. \end{equation}
The matrices
$\tilde {\boldsymbol{V}}$
have dimension
$d N \times ( N_t - d )$
, which can be large. However, the much larger spatial dimension underscores the importance of the dimensionality reduction step, which makes the implementation of the time-lagged snapshots feasible.
Then, another step of SVD removes potential redundancies introduced in (2.9)
\begin{equation} \tilde {\boldsymbol{V}}_1^{N_t - d + 1} \simeq \tilde {\boldsymbol{U}}\! \tilde {\boldsymbol{\varSigma }} \tilde {\boldsymbol{T}}^{\top } = \tilde {\boldsymbol{U}} \bar {\boldsymbol{T}}_1^{N_t-d+1}. \end{equation}
The number of modes is again truncated following
$\varepsilon _1$
, which determines the new spatial complexity
$N^{\prime } \leqslant N$
– if two thresholds are employed,
$\varepsilon _{1t}$
is used instead. Lastly, pre-multiplying (2.10) by
$\tilde {\boldsymbol{U}}$
and invoking (2.11) yields
\begin{equation} \bar {\boldsymbol{T}}_2^{N_t-d+1} = \bar {\boldsymbol{R}} \bar {\boldsymbol{T}}_1^{N_t-d}, \end{equation}
where
$\bar {\boldsymbol{R}} \simeq \tilde {\boldsymbol{U}}^{\top } \tilde {\boldsymbol{R}} \tilde {\boldsymbol{U}}$
has dimension
$N^{\prime } \times N^{\prime }$
.
2.2.3. Step 3: computation of the DMD modes
The DMD routine is applied to (2.12). The matrix
$\bar {\boldsymbol{T}}_1^{N_t-d}$
is decomposed using SVD (no truncation needed)
\begin{equation} \bar {\boldsymbol{T}}_1^{N_t-d} = \boldsymbol{U}\! \varSigma \boldsymbol{V}^{\top }. \end{equation}
Next,
$\bar {\boldsymbol{R}}$
is approximated such that
\begin{equation} \bar {\boldsymbol{R}} \simeq \bar {\boldsymbol{T}}_2^{N_t-d+1} \boldsymbol{V}\! \boldsymbol{\varSigma } \boldsymbol{U}^{\top }. \end{equation}
The spectral properties of
$\bar {\boldsymbol{R}}$
, characterised by its
$N^{\prime }$
eigenvalues
$\mu _m$
and eigenvectors
$\bar {\boldsymbol{q}}_m$
, provide information about the dynamics of the reduced snapshots, that can be reconstructed as
\begin{equation} \hat {\boldsymbol{v}}_k \simeq \sum _{m=1}^M a_m \boldsymbol{q_m} \mu _m^{k-1}, \quad k = 1, \ldots , N_t, \end{equation}
with the frequencies
$\omega _m$
and growth rates
$\delta _m$
in (2.6) given by
\begin{equation} \delta _m + i \omega _m = \textrm {log}\left ( \mu _m \right ) / \Delta t, \end{equation}
with
$\Delta t$
the (constant) sampling rate between temporal snapshots. The amplitudes of the DMD modes are computed using least-squares fitting in (2.15) (Chen, Tu & Rowley Reference Chen, Tu and Rowley2012). Finally, (2.6) is recovered from multiplying (2.15) by
$\boldsymbol{U}$
and invoking (2.8).
A second threshold
$\varepsilon _2$
defines the number of modes
$M$
, or spectral complexity, in (2.6), such that
\begin{equation} a_{M+1} / a_1 \leqslant \varepsilon _2. \end{equation}
By discarding low-amplitude modes, we obtain a sparse representation of the flow dynamics that only contains the most relevant structures. This indeed reduces the accuracy by neglecting small-scale features, and also prevents overfitting noise and spurious artefacts, thereby improving the generalisability of the reconstruction.
2.3. Spatio-temporal Koopman decomposition
The interpretability of the modes computed with HODMD can be improved by further decomposing them in space. This is the basis for the spatio-temporal Koopman decomposition (STKD) (Le Clainche & Vega Reference Le Clainche and Vega2018), that applies HODMD sequentially in time and space to approximate nonlinear data as a linear superposition of periodic or quasi-periodic standing and/or travelling waves.
To do so, HODMD is applied individually to each DMD mode in (2.6). For example, applying it in the spanwise direction yields
\begin{equation} \boldsymbol{u}_m \left ( x, y, z_r \right ) \simeq \sum _{n = 1}^{N} a_n \hat {\boldsymbol{u}}_{mn} \left ( x, y \right ) e^{\left ( \nu _{mn} i + \alpha _{mn} \right ) z_r}, \quad r = 1, \ldots , N_z, \end{equation}
with
$\kappa _{mn}$
and
$\alpha _{mn}$
the spanwise wavenumber and the growth rate, respectively, associated with each spanwise-temporal mode.
Finally, combining (2.18) with (2.6) gives
\begin{equation} \boldsymbol{v} ( x, y, z_r, t_k ) \simeq \sum _{m, n = 1}^{M, N} a_{mn} \hat {\boldsymbol{u}}_{mn} ( x, y ) e^{( \delta _m + \omega _m i ) t_k + ( \nu _{mn} + \alpha _{mn} i ) z_r}, k = 1, \ldots , N_t,\ r = 1, \ldots , N_z. \end{equation}
As a result, the data are now described as a superposition of travelling waves with spanwise-temporal amplitude
$a_{mn} = a_m a_n$
, whose phase velocity is defined as
$c_{mn} = \omega _m / \alpha _{mn}$
.
3. Global analysis
3.1. Data visualisation
Instantaneous fields of the streamwise fluctuating velocity
$u^{\prime }$
for the Newtonian (a) and the viscoelastic (b) turbulent planar jets. Panels a and e show an
$xy$
-plane extracted at
$z = 0$
. Insets provide a closer view of the near field up to
$x/h \approx 20$
. The dashed coloured lines indicate the
$xz$
-planes shown in panels b, c and d for the Newtonian and f, g and h for the viscoelastic jets.
Let us start by visually comparing the flow of the Newtonian and viscoelastic jets. Figure 2 shows an instantaneous field of the streamwise velocity fluctuation
$u^{\prime }$
with values comparable in the two cases despite the large difference in Reynolds number. Panels a and e are side views of the jets at
$z = 0$
, and we can observe that the far-field region is dominated by large-size structures in both jets. Here, while small-scale turbulence is also visible in the Newtonian jet, structures are smoother in the viscoelastic jet, with little evidence of small-scale fluctuations. This is consistent with the turbulent kinetic energy decaying significantly faster across scales in elastic turbulence than in Newtonian turbulence, with a signature energy spectrum decay in wavenumber of
$k^{-4}$
compared with the Kolmogorov one
$k^{-5/3}$
(Singh et al. Reference Singh, Perlekar, Mitra and Rosti2024). It is worth mentioning that, despite the exponent
$-4$
also being reported in viscoelastic planar jets at low Reynolds numbers (Rota et al. Reference Rota, Singh, Chiarini, Amor, Soligo, Mitra and Rosti2026), the energy spectrum in frequency decays as
$\omega ^{-3}$
(Yamani et al. Reference Yamani, Keshavarz, Raj, Zaki, McKinley and Bischofberger2021, Reference Yamani, Raj, Zaki, McKinley and Bischofberger2023; Soligo & Rosti Reference Soligo and Rosti2023); the difference between scaling exponents is a consequence of Taylor’s hypothesis not holding in elastic turbulence (Rota et al. Reference Rota, Singh, Chiarini, Amor, Soligo, Mitra and Rosti2026). The more dominant presence of large-size structures in the viscoelastic jet also influences the near field. Since the disturbance generated by the stronger coherent fluctuations can propagate upstream along the jet edges without being disrupted by the small-scale features, they reach even the proximity of the inlet, strongly affecting the onset of turbulence. This is not visible in the Newtonian case, where the small-scale fluctuations rapidly destroy their coherency and prevent large-scale noise propagating upstream.
Studying the near-field structures in more detail reveals the different nature of the flow instability (see insets in panels a and e). In the Newtonian planar jet, the flow is laminar with zero fluctuations close to the inlet, which start amplifying from around
$x/h \approx 5$
. Here, alternating positive–negative
$u^{\prime }$
lobes align along the upper and lower shear layers, consistently with Kelvin–Helmholtz rollers whose symmetry about the centreline matches the varicose (symmetric) instability, that dominates at high Reynolds numbers (Antonia et al. Reference Antonia, Browne, Rajagopalan and Chambers1983; Thomas & Goldschmidt Reference Thomas and Goldschmidt1986; Deo, Mi & Nathan Reference Deo, Mi and Nathan2008; Soligo et al. Reference Soligo, Chiarini and Rosti2025). The viscoelastic jet shows a notably different scenario. Streamwise velocity fluctuations rapidly appear in the form of streamwise-elongated structures along the shear layers of the potential core. This feature differs from the classical Kelvin–Helmholtz rollers, and indicates an alternative amplification pathway for the fluctuations, in which elasticity triggers an earlier transition – at lower Reynolds number and closer to the inlet – than in the Newtonian case at much higher Reynolds. The flow in both cases becomes unstable at larger distances from the inlet compared with some other works found in the literature (see for instance Da Silva & Métais (Reference Da Silva and Métais2002), Stanley et al. (Reference Stanley, Sarkar and Mellado2002), Guimarães et al. (Reference Guimarães, Pimentel, Pinho and da Silva2020)), since we do not add any disturbance at the inlet to induce an early transition.
Next, we focus on top views of the jet, with streamwise–spanwise planes extracted at three heights from the centreline, corresponding in the near field to the jet core (panels b and f), above the shear layer (panels c and g) and in the outer flow (panels d and h). In the near field of the Newtonian jet, spanwise-homogeneous Kelvin–Helmholtz rollers populate the shear layer, which destabilise and break-up moving downstream. In the far field, streamwise-elongated structures appear (see for instance panel b), reminiscent of the streaky structures reported by Nogueira et al. (Reference Nogueira, Cavalieri, Jordan and Jaunet2019) in round jets. These structures grow in size downstream as smaller structures merge. In the viscoelastic jet, roller-like structures are also present within the jet core (panel f), but their amplitude is much smaller than in the Newtonian case. In contrast, short and thin streaky structures populate the near-field region along the shear layer (panel g). They remain coherent at the potential core before they get disrupted by the large-scale fluctuations propagating upstream. These near-field streaks point toward a different, elasticity-driven pathway to turbulence in viscoelastic turbulent planar jets at low Reynolds numbers. Further downstream, the field is dominated by large-size, spanwise-coherent structures, with the flow remaining predominantly homogeneous in the spanwise direction and lacking the streamwise-elongated structures observed in the core of the Newtonian jet.
In the following sections, we aim to characterise the different structures observed here qualitatively, and investigate their dynamics and interactions.
3.2. Flow structures
The HODMD spectra overlapped for multiple calibrations for the Newtonian (upper panel) and viscoelastic (middle panel) jets. Non-dimensional frequency or Strouhal number,
$\textit{St} = f h / U$
, is compared with the normalised amplitude,
$a / a_1$
, with
$a_1$
being the largest amplitude in each series. Shaded areas indicate robust modes: blue refers to low-frequency modes or streamwise-elongated structures, while red denotes high-frequency modes or spanwise-coherent structures. The thickness of each bar matches the maximum deviation of the frequency in each cluster. Robust modes from both jets are compared in the lower panel (marker outlines are colour coded similar to the shaded areas).
We begin the characterisation of the global coherent structures in the Newtonian and viscoelastic jets using sequential HODMD or STKD. As a first step, data are decomposed in time to compute the temporal coherent structures. The wide range of temporal (and spatial) scales in both jets complicates the identification of flow patterns, so HODMD must be carefully calibrated to ensure physically meaningful results. Here, our interest is finding the large-amplitude modes associated with the dominant large-scale structures in the flow, which must be distinguished from the fairly large number of frequencies present in the data. To this end, we apply HODMD recursively with different combinations of
$\varepsilon _1$
,
$\varepsilon _2$
and
$d$
. Modes that are robust across calibrations, i.e. their frequency appears consistently regardless of parameter choice, are considered physical. In contrast, spurious modes arising from the mathematical decomposition tend to be scattered throughout the spectrum without clustering at particular frequencies.
This step is computationally expensive due to the recursive application of HODMD, so data are pre-processed to reduce their size. First, data are cropped in the streamwise and jet-normal directions, since the original simulations are performed in a much larger domain to avoid confinement effects. The cropped sub-domains have dimensions
$L_x = 140h$
and
$L_y = 60h$
; dimensions were chosen such that each box encloses a similar portion of the domain, which contains the jet and the fully developed flow. Second, the data are uniformly downsampled in all spatial directions by a factor of two. After this reduction, the Newtonian jet data have dimensions
$480 \times 207 \times 120$
, and the viscoelastic data
$337 \times 146 \times 64$
. Therefore, the HODMD is performed on the flow within the reduced box having the same lengths in the streamwise and jet-normal directions, and extending over the full span. The analysis is performed using a statistically stationary data set consisting of
$436$
snapshots for the Newtonian case and
$300$
for the viscoelastic case, that are sampled at the interval of
$\Delta t U / h = 2$
time units, that is sufficient for resolving the dynamics of the large-scale coherent structures in both jets.
The HODMD calibration is summarised in figure 3. In the viscoelastic jet, we set
$\varepsilon _1 = \varepsilon _2 = \varepsilon$
, whereas in the Newtonian one we split the first threshold in space and time,
$\varepsilon _{1s}$
and
$\varepsilon _{1t}$
, owing to the larger spatial complexity of the data. The calibration is then performed using
$\varepsilon = 6 \times 10^{-4}, 8 \times 10^{-4}$
and
$10^{-3}$
in the viscoelastic jet, and
$\varepsilon _{1t} = 10^{-4}, 2 \times 10^{-4}$
and
$3 \times 10^{-4}$
and
$\varepsilon _{1s} = 2 \times 10^{-3}$
in the Newtonian one, with
$d$
set to
$60, 70, 80, 90$
and
$100$
in both cases. We define two criteria for identifying robust modes: the first one over the frequency (the
$\omega$
-criterion) and the second one over the spatial shape of the DMD modes (the
$u$
-criterion). A mode with frequency
$\omega _m$
fulfils the
$\omega$
-criterion if
$|\omega _{m_i} - \omega _{m}| \leqslant \epsilon$
, with
$\epsilon = 5 \times 10^{-3}$
. If this criterion is fulfilled in
$75\,\%$
of the calibrations, modes are deemed common. Then, common modes are promoted to robust if their spatial shape
$\boldsymbol{u}_m$
is similar across calibrations, i.e. they fulfil the
$u$
-criterion. This is evaluated by computing the cosine similarity between DMD modes:
$\textrm {cos} ( \boldsymbol{u}_{m_i}, \boldsymbol{u}_{m_j} ) \geqslant \zeta$
, with
$\zeta = 0.8$
. We acknowledge robustness if
$\zeta$
is close to one, i.e. both modes are identical pairs; we request the
$u$
-criterion to be fulfilled in
$75\,\%$
of the calibrations. In doing this, we find twenty-one modes in the Newtonian jet and sixteen modes in the viscoelastic jet that are robust (the same number of modes were deemed common). Moreover, the amplitude across calibrations within each robust cluster is comparable (the deviation is of the order of
$10^{-2}$
at most), indicating that HODMD assigns a similar weight notwithstanding the choice of parameters.
In the following, we show the results for the set of parameters
$d = 80$
,
$\varepsilon _{1s} = 2 \times 10^{-3}$
and
$\varepsilon _{1t} = \varepsilon _2 = 2 \times 10^{-4}$
in the temporal analysis of the Newtonian jet, and
$d = 100$
and
$\varepsilon _1 = \varepsilon _2 = 6 \times 10^{-4}$
in the viscoelastic one. We chose the calibration that retrieved the largest number of robust modes, although a simpler criterion could be choosing an intermediate
$d$
and the lowest thresholds
$\varepsilon _1$
and
$\varepsilon _2$
. Figure 3(c) shows that the robust modes are predominantly at low and intermediate frequencies, with only a few high-frequency modes present in the Newtonian jet. In both jets, the dominant part of the spectra is shifted toward lower frequencies (modes are computed globally, so the spectra reflect the dominance of the low-frequency structures due to the energy distribution). The Newtonian spectrum, however, exhibits greater complexity. Smaller turbulent scales are hindered in the viscoelastic jet due to the low Reynolds number, while higher frequencies are more important for reconstructing the turbulent dynamics in the high-Reynolds-number Newtonian jet. On the other hand, the spectrum of the viscoelastic jet is dominated by more pronounced peaks, which resemble those reported by Suresh et al. (Reference Suresh, Srinivasan, Sundararajan and Das2008) in the power spectra of low-Reynolds-number Newtonian jets. In their case, the peaks corresponded with sub-harmonics of the fundamental vortex formation frequency; here, these are roughly harmonics of the dominant frequency
$\textit{St} \simeq 0.008$
(see for instance
$\textit{St} \simeq 0.023$
and
$0.056$
).
Spatial structure of robust DMD modes in the Newtonian and viscoelastic jets. Three-dimensional iso-surfaces represent the normalised streamwise velocity for magnitudes
$+0.5$
(red) and
$-0.5$
(blue). of the near field up to
$x/h \approx 30$
. The yellow semi-transparent surfaces mark the average jet thickness. Labels are colour coded similar to the shaded areas in figure 3 (blue for low-frequency or streamwise-elongated modes, and red for high-frequency or spanwise-coherent modes).
We next visualise the three-dimensional spatial shape
$\boldsymbol{u}_m$
of a few relevant robust modes in the Newtonian and viscoelastic jets in figure 4. Flow structures are broadly similar across low to intermediate frequencies between the cases: low-frequency modes (blue labels) represent streamwise-elongated structures, while higher-frequency modes show spanwise coherency (red labels). In the Newtonian jet, the lowest frequency (panel a) indicates streamwise-elongated structures that emerge just downstream of the potential core. These structures move along the edges of the jet, computed as the distance from the centreline where the local streamwise velocity equals half the value at the centreline, and they grow with downstream distance. In the viscoelastic jet, the slowest mode (panel f) also shows streamwise-elongated structures, although they are more localised compared with the same structures in the Newtonian case, which extend globally and have a dominant presence in the far field. Consistent with this picture, higher low-frequency modes (panels b and g) also capture the evolution of the streamwise-elongated structures. In the viscoelastic jet, this mode additionally highlights the short streaks in the potential core, similar to those in figure 2(g).
Higher frequencies show spanwise-coherent structures. In all cases, flow structures are confined rather than expanding globally like the low-frequency ones, and they are located symmetrically above and below the centreline of the jet. The shape and frequency of the modes at intermediate values (panels c in the Newtonian jet and h and i in the viscoelastic one) match the description of Orr wave packets (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018). At an even higher frequency, structures are displaced to the proximity of the inlet. In the Newtonian jet, these have enhanced spanwise homogeneity, for instance at
$\textit{St} \simeq 0.144$
(panel d) structures are located at the core and edges of the jet, the latter resembling the Kelvin–Helmholtz rollers. The spatial arrangement of the mode, symmetric about the centreline, and their frequency, very close to that reported experimentally by Deo et al. (Reference Deo, Mi and Nathan2008), for high-Reynolds-number Newtonian jets, is consistent with the symmetric (varicose) mode, also highlighted by Soligo et al. (Reference Soligo, Chiarini and Rosti2025) for the same database. It is noteworthy that, despite the method being applied globally, HODMD is also able to resolve the smaller-scale structures near the inlet. The method is also able to reconstruct very-high-frequency modes (panel e), that arise from nonlinear interactions between high-frequency modes at the near-field region. Similar structures are also highlighted at
$\textit{St} \simeq 0.086$
in the viscoelastic jet (panel j). In this case, spanwise-homogeneous rollers appear along the jet edges at much smaller frequency compared with the Newtonian jet. Despite the roll-up of the shear layer being significantly different between the Newtonian and viscoelastic jets (see insets in panels
$a$
and
$e$
in figure 2), the shape and location of both structures at high frequency (spanwise homogeneous and along the edges) point to Kelvin–Helmholtz rollers, although induced by a mechanism different from the one for Newtonian jets (Guimarães et al. Reference Guimarães, Pinho and da Silva2023). The same mode in the viscoelastic jet also represents small-scale structures along the jet edges downstream of the potential core, indicating the breakdown of the upstream rollers.
3.3. Spatio-temporal structures
The temporal analysis revealed the presence of similar structures in both planar jets: streamwise elongated at low frequency, and spanwise coherent at higher frequency. To complete their description, we now decompose each temporal DMD mode in the spanwise direction. For this analysis, we employ
$120$
snapshots, equidistant along the spanwise direction (
$\Delta z / h \simeq 0.29$
) in the Newtonian planar jet, and
$64$
snapshots, also sampled uniformly in spanwise (
$\Delta z / h \simeq 0.42$
) in the viscoelastic one, with thresholds
$\varepsilon _1 = \varepsilon _2 = 10^{-3}$
and
$d = 1$
in both cases. The decomposition yields an expansion of spatial modes that depend on the streamwise and jet-normal directions, and their superposition reconstructs each temporal DMD mode in terms of steady and travelling waves.
Spatio-temporal spectra. Normalised spanwise wavenumber,
$\kappa h$
, compared with their normalised spatio-temporal amplitude,
$\hat {a} / \hat {a}_{11}$
, with
$\hat {a}_{11}$
the largest amplitude in each series. Panels a and c correspond to low-frequency modes, while high-frequency modes are shown in b and d. A power-law decay of the normalised amplitude is suggested in each subpanel for high wavenumbers.
We show the spatial spectra for several robust frequencies in figure 5. Each amplitude quantifies the relative contribution of the different spanwise harmonics associated with each temporal mode. The smallest wavenumber is set by the spanwise domain length,
$L_z$
, such that
$\kappa _{\textit{min}} = 2\pi /L_z$
. In most cases, the spatio-temporal amplitudes decrease monotonically with the wavenumber, with the lowest non-zero wavenumber (
$\kappa _1 = \kappa _{\textit{min}}$
) carrying the largest contribution. The only exceptions are
$\textit{St} \simeq 0.029, 0.082$
in the Newtonian jet, where the second harmonic (
$\kappa _2 = 2 \kappa _1 = 2 \kappa _{\textit{min}}$
) dominates. The decay of the normalized amplitudes is similar, notwithstanding the frequency, with distinction made for the very-high-frequency modes (
$\textit{St} \simeq 0.144$
and
$0.206$
in the Newtonian jet, and
$\textit{St} \simeq 0.080$
and
$0.086$
in the viscoelastic jet), whose normalised amplitudes decay at a smaller rate. We also noticed two different decay rates, whether wavenumbers are above or below a threshold:
$\kappa \approx 0.7$
at low frequency (a) and
$\kappa \approx 1$
at high frequency (b) in the Newtonian jet, and
$\kappa \approx 2$
at low frequency (c) and
$\kappa \approx 3$
at high frequency (d) in the viscoelastic jet. The steeper decay of the normalised amplitudes in the viscoelastic jet highlights the greater complexity of the flow in the Newtonian jet.
Another sign of complexity is the low-rank behaviour. A large separation between the amplitudes at
$\kappa = 0$
and
$\kappa \gt 0$
indicates that the spatial mode with the largest amplitude dominates at a given frequency. This separation is more pronounced in the viscoelastic jet by almost one order of magnitude. For instance, at low frequency (panel c) the amplitude of the modes with wavenumber
$\kappa _1$
is significantly lower compared with those from the Newtonian jet (panel a), thus indicating that spanwise-homogeneous structures play a more dominant role in the slow dynamics of the viscoelastic jet, whereas spanwise-inhomogeneous (streamwise-elongated) structures are more relevant in the Newtonian jet (and subdominant in the viscoelastic one). At higher frequency, the amplitudes of the non-zero wavenumbers increase, and the low-rank behaviour reduces. However, at
$\textit{St} \simeq 0.144$
in the Newtonian jet (panel b) and at
$\textit{St} \simeq 0.080$
and
$0.086$
in the viscoelastic jet (panel d), the low-rank behaviour is the most pronounced (the amplitude of
$\kappa = 0$
is much larger than
$\kappa _1$
), hence the dominance of spanwise-homogeneous structures related to the Kelvin–Helmholtz instability at very high frequency (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020). The same low-rank behaviour is also observed at
$\textit{St} \simeq 0.023$
for the viscoelastic jet (panel d), although related in this case to homogeneous Orr wave packets that dominate at low frequency (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018).
Spatial structure of two low-frequency spatio-temporal modes. The upper panel shows
$xy$
-planes at
$z = 0$
for
$\kappa = 0$
, and the lower panel the three-dimensional iso-surfaces for
$\kappa _1$
(left column) and
$\kappa _2$
(right column). Both cases show the normalised streamwise velocity, with iso-surfaces indicating regions of magnitude
$+0.5$
(red) and
$-0.5$
(blue). The yellow lines and translucent surfaces mark the average jet thickness.
We complement the spatio-temporal spectra with the representation of the spatial shape of the spatio-temporal modes. We first show in figure 6 the reconstruction of the leading spatio-temporal modes for
$\kappa = 0$
, and
$\kappa _1$
and
$\kappa _2$
for two low-frequency modes. The
$\kappa = 0$
mode is spanwise homogeneous and therefore shown as a two-dimensional slice, while non-zero wavenumbers are represented using three-dimensional iso-surfaces. Note that the spatial modes are multiplied by
$e^{i \omega t} = \textrm {cos} (\omega t ) + i\textrm {sin} (\omega t )$
, i.e. they are complex, so they change their phase with time
$t$
. Therefore, each subpanel shows an instantaneous realisation of each mode, that for instance would be opposite if the phase of the mode is adjusted by
$\pi$
, i.e. positive iso-surfaces would become negative, and vice versa. At first, the modes with
$\kappa = 0$
highlight the dynamics of the fluid column in the Newtonian jet (orange panel). The magnitude of the modes is the largest at the potential core, where the flow is essentially two-dimensional. Downstream of the potential core, where the flow is turbulent, the magnitude of the mode rapidly decays and becomes nearly zero elsewhere. At the same location, the shape of the mode with
$\textit{St} \simeq 0.003$
in the viscoelastic jet (purple panel) is reminiscent of flapping oscillations, that are driven by large-scale coherent structures, that are represented at
$\textit{St} \simeq 0.008$
, indicated by the negative correlation between both halves of the jet from
$x/h \approx 50$
(Goldschmidt & Bradshaw Reference Goldschmidt and Bradshaw1973; Gortari & Goldschmidt Reference Gortari and Goldschmidt1981). In Newtonian jets, the antisymmetric mode dominates the dynamics at the far field from across Reynolds numbers (Antonia et al. Reference Antonia, Browne, Rajagopalan and Chambers1983; Thomas & Goldschmidt Reference Thomas and Goldschmidt1986). The same mode is more dominant at low Reynolds numbers (Deo et al. Reference Deo, Mi and Nathan2008; Suresh et al. Reference Suresh, Srinivasan, Sundararajan and Das2008; Soligo et al. Reference Soligo, Chiarini and Rosti2025), where the presence of larger structures causes a more vigorous mixing. Similar large-scale, coherent structures were also observed in high-Reynolds-number viscoelastic jets, promoting a more effusive large-scale stirring (Guimãraes, Pinho & Da Silva Reference Guimãraes, Pinho and Da Silva2025). This is also in good agreement with the faster spreading rate and centreline velocity decay reported by Soligo & Rosti (Reference Soligo and Rosti2023) in low-Reynolds-number viscoelastic jets when compared with Newtonian turbulent jets.
On the other hand, the modes with non-zero wavenumbers
$\kappa _1$
and
$\kappa _2$
show three-dimensional streamwise-elongated structures, that are similar to the streaky structures reported by Nogueira et al. (Reference Nogueira, Cavalieri, Jordan and Jaunet2019) and Pickering et al. (Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020) in Newtonian turbulent round jets. In the Newtonian jet (orange panel), large-scale streamwise-elongated structures have a global presence: they arrange downstream of the potential core along the jet edges, although they have great penetration within the jet core. Structures associated with
$\kappa _2$
are located closer to the potential core compared with those for
$\kappa _1$
. The modes in the viscoelastic jet (purple panel) have a spatial distribution similar to those of the Newtonian jet – after the potential core and along the jet edges – but they overall disappear after
$x/h \approx 100$
, they are less massive and they have little penetration within the jet core. This description correlates with the spatio-temporal spectra. The amplitude of the modes with
$\kappa _1$
and
$\kappa _2$
in the Newtonian jet is greater compared with those from the viscoelastic case, hence streamwise-elongated structures are more regular and persistent in the former. On the contrary, the streamwise-elongated structures in the viscoelastic jet are disrupted by the large-scale structures at the far field, thus breaking them and homogenising the flow.
Reconstruction of the near-field streaks in the viscoelastic jet. Real (a–d) and imaginary (panels e–h) components of the spatio-temporal modes with wavenumber
$\kappa = 12 \kappa _{1}$
. Three-dimensional iso-surfaces represent the normalised streamwise velocity of magnitude
$+0.5$
(red) and
$-0.5$
(blue). The black translucent dashed–dotted line indicates the centreline of the jet, and the red circle the coordinate
$x/h = 20$
.
While higher wavenumbers indicate smaller structures, they are significantly different between jets. At low frequency, these relate exclusively to the near-field streaks in the viscoelastic case, represented in figure 7. The smallest wavenumber associated with them is
$\kappa _{12} = 12 \kappa _{1}$
(wavenumbers greater than
$\kappa _2$
but smaller than
$\kappa _{12}$
correspond with streamwise-elongated structures located after the potential core), while its harmonics capture finer streaks in the same region. Across all low frequencies, this mode shows groups of streaks that are at most
$6h$
long, and they extend up to
$x/h \approx 20$
from the inlet along the upper and lower shear layers. These structures are well resolved (approximately three grid points thick), and they resemble the streaky pattern at the potential core in figure 2(g), with the mode capturing twenty-four alternating high- and low-speed streaks organised along the spanwise axis, consistent with the flow visualisation. The real and imaginary parts of the mode provide insight into the temporal dynamics of the streaks. Because the DMD modes are complex conjugates, their real and imaginary components represent the same structure shifted by half a period. Here, the two are not identical, suggesting that the flow reorganises within half a period. The streaks in the potential core therefore do not behave as standing structures, but they travel downstream and are disrupted locally in the streamwise direction.
To summarise, the spatio-temporal analysis of the low-frequency modes reveals large-scale, streamwise-elongated structures in the far field of both jets, and smaller-scale streaks in the near field of the viscoelastic jet. While the description of the far field is reasonably similar between cases, the presence of the near-field streaks influences the dynamics at the potential core in the viscoelastic jet. Here, high- and low-speed streaks induce regions of localised high shear and high elongation among approaching streaks, that stretch the polymer in the streamwise direction and ultimately induce the transition to turbulence due to elastic instabilities. This purely elastic pathway allows the transition to turbulence in viscoelastic planar jets at low Reynolds number and high polymer elasticity, and it differs from the stabilising effect at higher Reynolds number (Rallison & Hinch Reference Rallison and Hinch1995; Ray & Zaki Reference Ray and Zaki2015; Yamani et al. Reference Yamani, Keshavarz, Raj, Zaki, McKinley and Bischofberger2021, Reference Yamani, Raj, Zaki, McKinley and Bischofberger2023).
We conclude by characterising in figure 8 the spatial shape of four high-frequency modes in the Newtonian and viscoelastic jets. In the Newtonian jet (orange panel), the streamwise velocity field takes the form of wave packets at
$\kappa = 0$
. However, their placement changes from
$\textit{St} \geqslant 0.144$
: for lower
$\textit{St}$
, wave packets are placed downstream of the potential core at the upper and lower halves of the jet, while for higher
$\textit{St}$
they are at the potential core along the shear layers. Similarly, Pickering et al. (Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020) reported that the dominant mechanism at zero wavenumber in high-Reynolds-number Newtonian round jets changes from Orr to Kelvin–Helmholtz at
$\textit{St} \approx 0.3$
(based on the diameter of the nozzle), that is independent of the Reynolds number; for non-zero wavenumbers, the dominant response of the flow is determined by the Kelvin–Helmholtz mechanism. We do not find similitude with the Kelvin–Helmholtz wave packets for
$\kappa _1$
at
$\textit{St} \simeq 0.009$
and
$0.030$
but at
$\textit{St} \simeq 0.144$
and
$0.206$
– that are consistent with the description of Kelvin–Helmholtz wave packets – spanwise homogeneous and aligned along the shear layer.
In the viscoelastic jet (purple panel), the modes with
$\kappa = 0$
indicate Orr-type wave packets, with exception made at
$\textit{St} \simeq 0.086$
, that shows primarily the motion of the fluid column at the near field. On the other hand, the shape of the modes with wavenumber
$\kappa _1$
changes significantly with frequency: at
$\textit{St} \simeq 0.023$
, structures are rather streamwise coherent (Pickering et al. (Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020) also related non-zero wavenumbers at low frequency to modes of streaky structures) while the Orr wave packets are recovered at
$\textit{St} \simeq 0.034$
and
$0.056$
. Conversely,
$\textit{St} \simeq 0.086$
shows structures scattered throughout the domain, although those at the near field, located at the shear layers, resemble the Kelvin–Helmholtz wave packets.
4. Local analysis of the viscoelastic planar jet
In the previous section, we showed the existence of the near-field streaks, that potentially interact with the flow instability of the jet at the potential core. The previous analysis was global and thus it is difficult to fully assess the role of streaks in the near-field region. To address this, we now apply HODMD to the viscoelastic case in smaller three-dimensional boxes that enclose the flow near the inlet. Three boxes are considered with increasing streamwise length:
$l_x / h = 20, 30$
and
$40$
. All boxes are centred on the jet centreline, have a jet-normal length of
$l_y / h = 20$
and span the full width of the domain. In particular, the smallest box covers the potential core, while longer ones extend partially into the turbulent region.
Spatial structure of four high-frequency spatio-temporal modes. For each jet, the left column shows
$xy$
-planes at
$z = 0$
for
$\kappa = 0$
, and the right column the three-dimensional iso-surfaces for
$\kappa _1$
. Both cases show the normalised streamwise velocity, with iso-surfaces indicating regions of magnitude
$+0.5$
(red) and
$-0.5$
(blue). Insets provide a closer view at the near field up to
$x \approx 30h$
. The yellow lines and translucent surfaces mark the average jet thickness.
Local HODMD spectrum for robust frequencies. Non-dimensional frequency,
$\textit{St}$
, compared with the normalised amplitude,
$a / a_1$
, with
$a_1$
the largest amplitude in each series. Black solid circles mark frequencies detected in all three boxes, while black dashed circles indicate those present in two boxes.
Figure 9 shows the spectrum of robust modes in the three boxes. The method was recalibrated for the local analysis:
$d$
is varied from
$60$
to
$120$
(larger
$d$
allows us to capture a more complex dynamics because of the higher dimension of the delay-embedding space), while the threshold
$\varepsilon$
is kept the same as the global analysis, with the addition of
$\varepsilon = 4 \times 10^{-4}$
to spot smaller-amplitude modes. In the following, we show the results for
$d = 70$
in the smallest box and
$d = 80$
in the remaining two, while
$\varepsilon _1$
and
$\varepsilon _2$
are set to
$6 \times 10^{-4}$
in all cases. The smallest box yields
$7$
robust modes, compared with
$16$
and
$14$
for the longer ones. The dominant mode also shifts to lower frequencies with distance –
$\textit{St} \simeq 0.086$
in the smallest box,
$\textit{St} \simeq 0.056$
in the middle-sized one and
$\textit{St} \simeq 0.034$
in the longest box – consistent with the growing influence of a slower, larger-scale dynamics farther from the potential core. Similarly, the cumulative amplitude of the high-frequency modes is higher in the smallest box since scales are faster and smaller, so high frequencies are more relevant. We should point out that some of these modes likely persist in the longer boxes, but their amplitude may fall below the threshold
$\varepsilon$
used in the analysis, and are thus not detected by the method. Moreover, we identify a mode with very small frequency in the longer boxes, although this is not the case for the smallest one, where the dynamics associated with lower frequencies (streaks) is hindered in the smallest box because of the dominance of smaller-scale, higher-frequency modes. Although some of the high-frequency modes, in particular those in the smallest box, might be related to the jet instability, none of the extracted frequencies show temporal growth – the jets are stable in time – thus indicating that the instabilities are not absolute in nature.
Local spatio-temporal spectra. Normalised spanwise wavenumbers,
$\kappa h$
, compared with their normalised amplitude,
$\hat {a} / \hat {a}_{11}$
, from robust modes computed in all three boxes. Similar frequencies have the same colours among plots.
We extend the temporal analysis by decomposing the modes in the spanwise direction. Figure 10 shows the spatio-temporal spectra obtained in all three boxes. Overall, the amplitude of the non-zero wavenumbers increases with box length, highlighting the growing importance of three-dimensional structures in the turbulent region. The dominant wavenumber also shifts from
$\kappa _2 = 2 \kappa _{\textit{min}}$
in the smallest box to
$\kappa _1 = \kappa _{\textit{min}}$
in the longer ones, indicating that larger structures emerge downstream of the potential core. In the two longer boxes, the lowest-frequency mode has the smallest amplitude, so streamwise-elongated structures are sub-dominant in the near field. Moreover, the dominant wavenumber is
$\kappa _2$
, in good agreement with the global analysis, which suggested that smaller streamwise-elongated structures, and therefore higher wavenumbers, are located closer to the potential core. This trend is similar at higher frequency, with
$\kappa _1$
modes increasing their amplitude the most with box length. In particular,
$\kappa _1$
at
$\textit{St} \simeq 0.081$
and
$0.086$
grows significantly, eventually surpassing
$\textit{St} \simeq 0.056$
in the longest box. Finally, the amplitude of
$\textit{St} \simeq 0.121$
and
$0.164$
decreases with the box length (
$\textit{St} \simeq 0.121$
is not found outside the smallest box), where
$\kappa _2$
for
$\textit{St} \simeq 0.121$
dominates the spectrum in the smallest box.
Reconstruction of the bulk flow in the near field. Two-dimensional
$xy$
-planes at
$z = 0$
of the normalised streamwise velocity component for the modes with
$\kappa = 0$
.
Next, we visualise the spatio-temporal modes in all three boxes. We first show the modes with
$\kappa = 0$
in figure 11. In the global analysis, these modes were linked to two-dimensional wave packets. Here, we find qualitatively similar structures in the longest box, while the smaller boxes show the motion of the bulk flow. Starting with the longest box,
$\kappa = 0$
at the lowest frequency
$\textit{St} \simeq 0.007$
matches the mean flow, and it shows the dissipation of the jet downstream of the potential core as it spreads. At higher frequency, wave packet-like structures appear. They are located in the upper and lower halves of the jet for
$\textit{St} \simeq 0.034$
and
$0.056$
, while they align with the centreline for
$\textit{St} \simeq 0.086$
, while the flow structures at
$\textit{St} \simeq 0.081$
resemble a combination of the mentioned three modes. Looking to the smaller boxes for the same frequencies, we can connect these patterns to the antisymmetric or sinuous mode (
$\textit{St} \simeq 0.034, 0.056$
) and the symmetric or varicose mode (
$\textit{St} \simeq 0.086$
) of the bulk flow. In the smallest boxes, the modes clearly show the sinuous mode for
$\textit{St} \simeq 0.056$
and the varicose one for
$\textit{St} \simeq 0.086$
and
$0.104$
, that are visible from
$x \approx 10h$
. Even though both modes are obtained in the viscoelastic jet, the sinuous mode is amplified the most downstream, as indicated by the growth of the magnitude of the mode, that also correlates with the growth of its temporal amplitude in figure 9. Similarly, the mode at
$\textit{St} \simeq 0.086$
experiences a similar growth in space, although it is confined to
$10h \lessapprox x \lessapprox 30h$
. Interestingly, for
$\textit{St} \simeq 0.121$
, we observe fluctuations mainly in the shear layer rather than in the jet column, that destabilise the flow downstream. Yamani et al. (Reference Yamani, Raj, Zaki, McKinley and Bischofberger2023) reported a similar transition induced by a shear layer instability in elasto-inertial turbulent planar jets. At high Reynolds number, this instability arises from the interplay between inertia and elasticity, driven by the crowding and dilation of streamlines within the edge region of the jet (Rallison & Hinch Reference Rallison and Hinch1995). In elastic turbulence, inertia is negligible, and high shear originates from the near-field streaks instead, as shown in the previous section.
Reconstruction of the wave packet modes at the near field. Three-dimensional iso-surfaces of the modes with
$\kappa _1$
(a–c) and
$\kappa _2$
(d–f) of the normalised streamwise velocity for values of
$+0.5$
(red) and
$-0.5$
(blue). The background corresponds to
$xy$
-planes at
$z = 0$
of the mode with
$\kappa = 0$
for the same frequency. The yellow translucent surfaces mark the average jet thickness.
The results for
$\kappa = 0$
reveal that the near-field streaks not only modify the bulk flow, but they also interact with the coherent structures in the same region. This interaction ultimately drives their breakdown, as illustrated in figure 7. To study this in more detail, we show in figure 12 the three-dimensional modes with
$\kappa _1$
and
$\kappa _2$
for
$\textit{St} \simeq 0.056$
,
$0.086$
and
$0.121$
, that have the largest amplitudes. In all cases, the modes resemble wave packet structures. At
$\textit{St} \simeq 0.056$
, these structures are in the shear layer, whereas at
$\textit{St} \simeq 0.086$
and
$0.121$
they appear in both the shear layer and jet column for different wavenumbers. We therefore distinguish two types of modes: a shear layer mode at higher wavenumber and a dual mode (shear layer and jet column) mode at lower wavenumber. The spatio-temporal analysis separates their contribution and quantifies their importance based on the relative spatio-temporal amplitude: the dual mode is more relevant at
$\textit{St} \simeq 0.086$
, while the shear layer mode dominates at
$\textit{St} \simeq 0.121$
. Moreover, the shear layer mode develops earlier in the potential core than the dual mode – more specifically the structures in the jet column – and its cumulative amplitude across frequencies exceeds that of the dual mode. Taken together, these findings indicate that the shear layer dynamics plays a more significant role at the potential core. This classification aligns with theoretical and experimental results. The shear layer mode generally appears at higher wavenumbers than the jet column mode – the dual mode in this case – and they merge at sufficiently high elasticity (Rallison & Hinch Reference Rallison and Hinch1995; Yamani et al. Reference Yamani, Raj, Zaki, McKinley and Bischofberger2023). Finally, both modes, especially the shear layer one, interact with the near-field streaks. The interaction is the most relevant at
$\textit{St} \simeq 0.086$
and
$0.121$
, owing to the earlier onset of the structures. The near-field streaks also interact with the mode at
$\textit{St} \simeq 0.056$
, although such interaction occurs near the end of the potential core at
$x/h \approx 20$
between the wave packets and the remnants of the near-field streaks that broke upstream. This interplay between streaks and wave packets constitutes the mechanism responsible for streak breakdown in the near field, and it completes the effect of the streaks – modification of bulk flow and interaction with flow instability – at the potential core. This interplay is elastic in nature, where both near-field streaks and wave packets are driven by the elasticity of the flow.
5. Polymer dynamics
We conclude our analysis by examining the trace of the conformation tensor to highlight similarities and differences between coherent structures in the velocity and polymer fields of the viscoelastic planar jet. The spatio-temporal analysis is applied globally to the three-dimensional trace of the conformation tensor. The calibration from STKD is analogous to that used for the velocity field and, in particular, we show the results for the set of parameters
$d = 80$
and
$\varepsilon _1 = \varepsilon _2 = 6 \times 10^{-4}$
.
Temporal (a) and spatio-temporal (b) spectra of the trace of conformation tensor. Panel a shows the non-dimensional temporal frequency,
$\textit{St}$
, compared with the normalised amplitude,
$a / a_1$
, with
$a_1$
the largest amplitude in each series. The robust modes from the polymer field (light purple) are compared with those from the velocity field (purple) reported in figure 3. Red outlines indicate modes of the polymer that are promoted to the spatio-temporal analysis; panel b shows the normalised spanwise wavenumber,
$\kappa h$
, compared with the normalised spatio-temporal amplitude,
$\hat {a} / \hat {a}_{11}$
, with
$\hat {a}_{11}$
the largest amplitude for each robust mode. A power-law decay of the normalised amplitude is suggested for high wavenumbers.
Figure 13 summarises the spatio-temporal spectra. Panel a shows the temporal frequencies and amplitudes of the robust modes, that are compared with those from the velocity field. Panel b presents the corresponding spanwise spectra for the subset of robust modes highlighted with red outlines in panel a. The velocity and polymer fields share a similar temporal dynamics: several robust frequencies identified in the analysis of the velocity field also appear in the corresponding for the polymer field. However, the polymer strongly amplifies these modes, indicated by their higher amplitude, that are roughly an order of magnitude larger. This amplification is particularly evident at higher frequencies, where the dominant frequency shifts to
$\textit{St} \simeq 0.026$
, suggesting that faster dynamics is more relevant for reconstructing the polymer field. The spanwise spectrum also resembles that of the velocity field. In both cases, amplitudes decay monotonically with increasing wavenumber, indicating that smaller scales contribute progressively less to the reconstruction of the temporal dynamics. For the polymer, however, this decay is less pronounced, as indicated by smaller exponent of the power-law decay of the normalised amplitude (
$-1$
compared with
$-1.5$
in the velocity field), reflecting a more complex structure. Moreover, the spatio-temporal amplitudes of the non-zero wavenumbers at low frequency are higher, indicating that the streaky structure of the polymer field is more regular compared with the velocity field.
Spatial structures of two low-frequency spatio-temporal modes from the polymer field. Panels a and b show the modes with
$\kappa _1$
, and panels c and d with
$\kappa _2$
. Three-dimensional iso-surfaces are plotted for magnitudes
$+0.5$
(red) and
$-0.5$
(blue). The yellow translucent surfaces mark the average jet thickness.
We now compare the spatio-temporal structures of the polymer with those from the velocity field. First, we show in figure 14 the two leading non-zero wavenumbers,
$\kappa _1$
and
$\kappa _2$
, at low frequency. The first distinctive structures of the polymer field are polymer filaments stretched in the streamwise direction. These filaments are connected to the same low frequencies of the velocity streaks, exhibit a qualitatively similar spatial organisation and populate the same region of the flow. At higher wavenumber (panels c and d), smaller filaments are located more upstream, similar than the velocity streaks, although they populate a region closer to the inlet (panel d). The size of these structures is larger compared with the velocity streaks at the same location. At the potential core, the strong shear generated from high- and low-speed streaks stretches the polymer; the only structures capable of stretching the polymer are indeed the near-field streaks. This confirms that polymer stretching occurs precisely in this region, supporting the view that elasticity plays a decisive role in triggering the transition to elastic turbulence.
Spatial structures of three high-frequency spatio-temporal modes from the polymer field. Two-dimensional
$xy-$
planes at
$z = 0$
of the normalised streamwise velocity component for the modes with
$\kappa = 0$
. The yellow lines mark the average jet thickness.
Lastly, the second distinctive structures of the polymer field are the centre-mode structures. These are shown in figure 15. In this case, polymer sheets first appear in both shear layers, are stretched by the faster flow in the jet core and then merge at the centre further downstream. At
$\textit{St} \simeq 0.056$
and
$0.081$
, stretching occurs immediately after the potential core, and the resulting structures are smaller than those observed at
$\textit{St} \simeq 0.026$
. The organisation also changes with frequency. While the structure is antisymmetric with respect to the centreline at
$\textit{St} \simeq 0.026$
, it becomes symmetric from
$\textit{St} \simeq 0.056$
. Similar highly extended sheets of polymer stress have been reported in two-dimensional elasto-inertial channel flows (Dubief, Terrapon & Hof Reference Dubief, Terrapon and Hof2023). These are connected to the centre-mode instability of viscoelastic pipe and channel flows at high Reynolds number (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021), that persist at vanishing inertia (Buza et al. Reference Buza, Beneitez, Page and Kerswell2022). While similar centre-mode structures have been observed in elastic turbulent channel flows (Morozov Reference Morozov2022; Rota et al. Reference Rota, Amor, Le Clainche and Rosti2024), they also appear in free-shear flow configurations (Berti et al. Reference Berti, Bistagnino, Boffeta, Celani and Musacchio2008; Lewy & Kerswell Reference Lewy and Kerswell2025). Most of these observations were done in a two-dimensional domain; in fact, the centre-mode structure of the viscoelastic planar jet is two-dimensional, where the superposition of
$\kappa = 0$
with the leading non-zero wavenumber and its harmonics modulates the two-dimensional structure in the three-dimensional space. While it is believed that elasto-inertial turbulence is essentially a two-dimensional phenomenon (Sid, Terrapon & Dubief Reference Sid, Terrapon and Dubief2018; Dubief et al. Reference Dubief, Terrapon and Hof2023), elastic turbulence in planar jets is necessarily three-dimensional: the near-field streaks, that are inhomogeneous structures in the spanwise direction, play a relevant role in the transition to elastic turbulence. These will be unavoidably missed if the domain is fully two-dimensional, remarking the three-dimensional nature of elastic turbulence.
6. Conclusions
We perform a data-driven modal analysis of Newtonian and viscoelastic turbulent planar jets using STKD, where we represent the flow dynamics as a superposition of modes related to the large-scale coherent structures. This sparse representation facilitates the interpretation of the complex dynamics in turbulent jet flows, allowing us to compare the most relevant structures between the high-Reynolds-number Newtonian and the low-Reynolds-number viscoelastic jets. We therefore compare a fully developed Newtonian turbulent jet with a fully developed elastic turbulent jet; in the latter case, turbulence is sustained solely by elasticity, owing to the low Reynolds number and high Deborah number, with the equivalent low-Reynolds-number Newtonian jet being laminar (Sato Reference Sato1960; Sato & Sakao Reference Sato and Sakao1964; Soligo & Rosti Reference Soligo and Rosti2023).
The global analysis of the velocity field reveals the presence of slow streaky structures and fast wave packets. While their overall characterisation is similar, there are subtle differences. Low-frequency streaky structures have a global presence in the Newtonian jet, particularly in the far field, where they also extend deep into the jet core. In contrast, they are less dominant in the viscoelastic jet, where the strong flapping motion at the wake disrupts the far-field streaks, rendering the flow more uniform and homogeneous. High-frequency modes show closer agreement between the two cases, with structures consistent with the description of Orr and Kelvin–Helmholtz wave packets. Moreover, the global analysis captures the characteristic instability in the Newtonian jet, where symmetries in the high-frequency modes align with the symmetric, varicose instability, known to dominate at high Reynolds numbers (Antonia et al. Reference Antonia, Browne, Rajagopalan and Chambers1983; Thomas & Goldschmidt Reference Thomas and Goldschmidt1986; Suresh et al. Reference Suresh, Srinivasan, Sundararajan and Das2008; Deo et al. Reference Deo, Mi and Nathan2008). In the viscoelastic jet, the wave packets more closely resemble the Orr wave packets, while evidence of Kelvin–Helmholtz wave packets, particularly at intermediate frequencies and non-zero wavenumbers, is more inconclusive.
The most striking difference between the two jets is the presence of streaks in the near field of the viscoelastic jet, that are absent in the Newtonian case and arise from flow elasticity. The near-field streaks are not standing structures: they travel and get disrupted by the flow instability of the jet. Local analysis of the potential core shows that the streaks interact with the high-frequency modes in the potential core. In particular, we identified two modes: a low-wavenumber dual mode with wave packets located in both the jet column and the shear layer, and a high-wavenumber shear layer mode. Similar modes have been reported by Yamani et al. (Reference Yamani, Raj, Zaki, McKinley and Bischofberger2023) in elasto-inertial turbulent planar jets using DMD. These modes are elastic in nature – inertial effects are negligible in the viscoelastic jet – where their interaction with the streaks induces the breakdown of the latter. At the same time, the streaks modify the bulk flow motion at the potential core. We observe both a symmetric and an antisymmetric dynamics: the symmetric motion dominates in the potential core, while the antisymmetric motion get amplified further downstream. The bulk flow is also perturbed at the shear layer, where disturbances emerge immediately after the inlet, grow downstream and ultimately destabilise the jet. These observations are consistent with the transition reported along the jet edges in viscoelastic jets at higher Reynolds numbers (Rallison & Hinch Reference Rallison and Hinch1995; Yamani et al. Reference Yamani, Keshavarz, Raj, Zaki, McKinley and Bischofberger2021, Reference Yamani, Raj, Zaki, McKinley and Bischofberger2023). Therefore, we link the near-field streaks to the onset of elastic turbulence in viscoelastic planar jets. We propose that high- and low-speed streaks generate regions of localised shear between adjacent streaks in the near field of the jet. The high shear stretches the polymer in the streamwise direction, the resulting elastic stresses are amplified downstream and they ultimately cause the transition to turbulence. Our findings thus point to the near-field streaks as a potential pathway to elastic turbulence in low-Reynolds-number viscoelastic planar jets.
Finally, we applied the same analysis to the polymer field. The temporal analysis showed that high-frequency modes are more important for reconstructing the polymer dynamics, which mainly reflects the amplified response to the velocity field. As a result, the polymer field exhibits similar structures: streamwise-oriented filaments of stretched polymer at low frequency and centre-mode structures at high frequency. The polymer filaments correspond to the same low frequencies as the velocity streaks, share a similar spatial organisation and occupy the same region of the flow. Notably, they also appear closer to the inlet at much smaller wavenumbers than the velocity streaks. This behaviour suggests that strong shear caused by high- and low-speed streaks stretches the polymer, reinforcing the idea that the polymer dynamics contributes actively to the transition. Together, velocity streaks and polymer filaments highlight that elastic turbulence is inherently three-dimensional in viscoelastic planar jets, similarly to what was observed recently in planar wall-bounded elastic turbulent flows (Lellep et al. Reference Lellep, Linkmann and Morozov2024; Rota et al. Reference Rota, Amor, Le Clainche and Rosti2024). In the planar jet, neglecting the inhomogeneous streaks and filaments would overlook key mechanisms driving the transition. On the other hand, we also identified high-frequency, centre-mode structures that resemble the arrowheads or narwhals reported in elasto-inertial and elastic turbulence in wall-bounded and free-shear flows (Berti et al. Reference Berti, Bistagnino, Boffeta, Celani and Musacchio2008; Morozov Reference Morozov2022; Dubief et al. Reference Dubief, Page, Kerswell, Terrapon and Steinberg2022; Rota et al. Reference Rota, Amor, Le Clainche and Rosti2024; Lewy & Kerswell Reference Lewy and Kerswell2025). Their presence confirms that the viscoelastic planar jet attains the elastic turbulent regime.
Acknowledgements
The authors acknowledge the Scientific Computing & Data Analysis section of the Core Facilities at OIST and HPCI for the provision of computational resources under Research Project grants hp220099, hp230018 and hp250035. Part of this work was carried out during the 2023 Madrid Turbulence Workshop, organised by Professor J. Jiménez and supported by ERC through Caust Grant ERC-AdG-101018287.
Funding
This research was supported by the Okinawa Institute of Science and Technology Graduate University (OIST) through subsidy funding to M.E.R. from the Cabinet Office, Government of Japan. M.E.R. also acknowledges support from the Japan Society for the Promotion of Science (JSPS) under grants 24K17210 and 24K00810. S.L.C. acknowledges support from grant PID2023-147790OB-I00 funded by MCIU/AEI/10.13039/501100011033/FEDER, UE.
Declaration of interests
The authors report no conflict of interest.
Open access
